{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sympy as sy\n",
    "M = sy.MatrixSymbol('M', 3, 3)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}a & 2 & 2 a\\\\2 & 1 & b\\\\b & a & 3\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([\n",
       "[a, 2, 2*a],\n",
       "[2, 1,   b],\n",
       "[b, a,   3]])"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from sympy import Symbol\n",
    "a = Symbol('a')\n",
    "b = Symbol('b')\n",
    "A = sy.Matrix(3,3,[a,2,b,2,1,a,2*a,b,3])\n",
    "B=A.transpose()\n",
    "\n",
    "B[1,1]\n",
    "B"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle a + 4$"
      ],
      "text/plain": [
       "a + 4"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sy.trace(B)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'unit_x' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[0;32m/tmp/ipykernel_9116/133189190.py\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m     40\u001b[0m \u001b[0mT_B_C\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mtransformation_matrix_qt\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'q_B_C_x'\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m'q_B_C_y'\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m'q_B_C_z'\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m'q_B_C_w'\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m't_B_C_x'\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m't_B_C_y'\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m't_B_C_z'\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     41\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 42\u001b[0;31m \u001b[0mA\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mT_W_B\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0mT_B_C\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0munit_x\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     43\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     44\u001b[0m \u001b[0;32mfor\u001b[0m \u001b[0mi\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mrange\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mA\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mshape\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mNameError\u001b[0m: name 'unit_x' is not defined"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "def transformation_matrix_qt(q,t):\n",
    "    qx=Symbol(q[0])\n",
    "    qy=Symbol(q[1])\n",
    "    qz=Symbol(q[2])\n",
    "    qw=Symbol(q[3])\n",
    "\n",
    "    px=Symbol(t[0])\n",
    "    py=Symbol(t[1])\n",
    "    pz=Symbol(t[2])\n",
    "\n",
    "    \n",
    "    \n",
    "    qx_2=Symbol(q[0]+\"*\"+q[0])\n",
    "    qy_2=Symbol(q[1]+\"*\"+q[1])\n",
    "    qz_2=Symbol(q[2]+\"*\"+q[2])\n",
    "    \n",
    "    T=sy.Matrix(4,4,[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0])\n",
    "     \n",
    "    T[0,0]=1-2*qy_2-2*qz_2\n",
    "    T[1,1]=1-2*qx_2-2*qz_2\n",
    "    T[2,2]=1-2*qx_2-2*qy_2\n",
    "\n",
    "    T[0,1]=2*(qx*qy+qw*qz)\n",
    "    T[2,0]=2*(qy*qz+qw*qx)\n",
    "    T[1,2]=2*(qx*qz+qw*qy)\n",
    "\n",
    "    T[0,2]=2*(qx*qz-qw*qy)\n",
    "    T[1,0]=2*(qx*qy-qw*qz)\n",
    "    T[2,1]=2*(qy*qz-qw*qx)\n",
    "\n",
    "    T[0,3]=px\n",
    "    T[1,3]=py\n",
    "    T[2,3]=pz\n",
    "    T[3,3]=1\n",
    "    \n",
    "    return T\n",
    "\n",
    "T_W_B=transformation_matrix_qt(['q_x','q_y','q_z','q_w'],['p_x','p_y','p_z'])\n",
    "T_B_C=transformation_matrix_qt(['q_B_C_x','q_B_C_y','q_B_C_z','q_B_C_w'],['t_B_C_x','t_B_C_y','t_B_C_z'])\n",
    "\n",
    "A=T_W_B*T_B_C*unit_x\n",
    "\n",
    "for i in range(A.shape[0]):\n",
    "    for j in range(A.shape[1]):\n",
    "        print(\"P[\"+str(i)+\",\"+str(j)+\"]=\\n\"+str(A[i,j]))\n",
    "        \n",
    "        print(\"   \")\n",
    "        print(\"  \")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "$\\left[\\begin{matrix}\\frac{\\left(p_{T x} - p_{x}\\right) \\cos{\\left(\\theta \\right)} + \\left(p_{T y} - p_{y}\\right) \\sin{\\left(\\theta \\right)}}{\\sqrt{\\left|{p_{T x} - p_{x}}\\right|^{2} + \\left|{p_{T y} - p_{y}}\\right|^{2} + \\left|{p_{T z} - p_{z}}\\right|^{2}}}\\end{matrix}\\right]$\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "\n",
    "\n",
    "\n",
    "def transformation_matrix_dh(theta,alpha,a_,d_):\n",
    "    st=Symbol(\"sin(\"+theta+\")\")\n",
    "    ct=Symbol(\"cos(\"+theta+\")\")\n",
    "    sa=Symbol(\"sin(\"+alpha+\")\")\n",
    "    ca=Symbol(\"cos(\"+alpha+\")\")\n",
    "    \n",
    "    st=sy.sin(theta)\n",
    "    ct=sy.cos(theta)\n",
    "    sa=sy.sin(alpha)\n",
    "    ca=sy.cos(alpha)\n",
    "    \n",
    "    a=Symbol(a_)\n",
    "    d=Symbol(d_)\n",
    "    \n",
    "    T=sy.Matrix(4,4,[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0])\n",
    "    \n",
    "    T[0,0]=ct\n",
    "    T[1,1]=ct*ca\n",
    "    T[2,2]=ca\n",
    "\n",
    "    T[0,1]=-st*ca\n",
    "    T[2,0]=0\n",
    "    T[1,2]=-ct*sa\n",
    "\n",
    "    T[0,2]=st*sa\n",
    "    T[1,0]=st\n",
    "    T[2,1]=sa\n",
    "\n",
    "    T[0,3]=a*ct\n",
    "    T[1,3]=a*st\n",
    "    T[2,3]=d\n",
    "    T[3,3]=1\n",
    "    \n",
    "    return T\n",
    "\n",
    "trans=transformation_matrix_dh(\"theta\",\"alpha\",\"a\",\"d\")\n",
    "    \n",
    "unit_x=sy.Matrix(4,1,[1,0,0,1])\n",
    "unit_y=sy.Matrix(4,1,[0,1,0,1])\n",
    "unit_z=sy.Matrix(4,1,[0,0,1,1])\n",
    "\n",
    "\n",
    "\n",
    "trans.subs('alpha',0)\n",
    "\n",
    "unitX=sy.Matrix(3,1,[1,0,0])\n",
    "\n",
    "px,py,pz=sy.symbols(['p_x','p_y','p_z'])\n",
    "ptx,pty,ptz=sy.symbols(['p_T_x','p_T_y','p_T_z'])\n",
    "\n",
    "pb=sy.Matrix(3,1,[px,py,pz])\n",
    "pt=sy.Matrix(3,1,[ptx,pty,ptz])\n",
    "\n",
    "vtb=pt-pb\n",
    "\n",
    "dist=vtb.norm()\n",
    "\n",
    "vc=trans[:3,:3]*(unitX)\n",
    "\n",
    "\n",
    "print(\"$\"+sy.latex(vc.transpose()*vtb/dist)+\"$\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$\\left[\\begin{matrix}\\frac{\\left(p_{T x} - p_{x}\\right) \\cos{\\left (\\theta \\right )} + \\left(p_{T y} - p_{y}\\right) \\sin{\\left (\\theta \\right )}}{\\sqrt{\\left|{p_{T x} - p_{x}}\\right|^{2} + \\left|{p_{T y} - p_{y}}\\right|^{2} + \\left|{p_{T z} - p_{z}}\\right|^{2}}}\\end{matrix}\\right]$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sympy as sy\n",
    "\n",
    "d=sy.symbols('d')\n",
    "\n",
    "def Roderiguez(v1,v2):\n",
    "    \n",
    "    v1n=v1.normalized()\n",
    "    v2n=v2.normalized()\n",
    "    \n",
    "    k=v1n.cross(v2n).normalized()\n",
    "    ct=(v1n.transpose()*v2n)[0]\n",
    "    \n",
    "    st=sy.sqrt(1-ct*ct)\n",
    "    kx=k[0]\n",
    "    ky=k[1]\n",
    "    kz=k[2]\n",
    "    \n",
    "    k_as=sy.Matrix(3,3,[0,-kz,ky,kz,0,-kx,-ky,kx,0])\n",
    "    \n",
    "    I=sy.Matrix(3,3,np.eye(3).flatten())\n",
    "    \n",
    "    R=I*ct+(1-ct)*k*(k.transpose())+st*k_as\n",
    "    \n",
    "    \n",
    "    return R\n",
    "    \n",
    "v1=vc.subs('theta',0)\n",
    "\n",
    "x=np.array([[-2.01,-1.001]])\n",
    "\n",
    "sub_key=['p_x','p_y','p_z','p_T_x','p_T_y','p_T_z']\n",
    "sub_val=[x[0,0],x[0,1],0,0,0,0]\n",
    "sub_dict=dict(zip(sub_key,sub_val))\n",
    "\n",
    "v2=vtb.subs(sub_dict)\n",
    "R=Roderiguez(v1,v2)\n",
    "\n",
    "s_x,s_y,s_z=sy.symbols(['sx','sy','sz'])\n",
    "S=sy.Matrix(3,3,np.diag([s_x,s_y,s_z]).flatten())\n",
    "\n",
    "#sub_key=['sx','sy','sz']\n",
    "#sub_val=[10,1,2]\n",
    "#sub_dict=dict(zip(sub_key,sub_val))\n",
    "\n",
    "P=(R*S*R.transpose())#.subs(sub_dict)\n",
    "\n",
    "#P_np=np.array([[P[0,0],P[0,1]],[P[1,0],P[1,1]]],dtype=np.double)\n",
    "#print(\"$\"+sy.latex(P[0,0])+\"$\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$sx \\left(\\frac{1.0 p_{T x} \\cos{\\left (\\theta \\right )}}{\\sqrt{\\left|{\\sin{\\left (\\theta \\right )}}\\right|^{2} + \\left|{\\cos{\\left (\\theta \\right )}}\\right|^{2}} \\sqrt{\\left|{p_{T x}}\\right|^{2} + \\left|{p_{T y}}\\right|^{2} + \\left|{p_{T z}}\\right|^{2}}} + \\frac{1.0 p_{T y} \\sin{\\left (\\theta \\right )}}{\\sqrt{\\left|{\\sin{\\left (\\theta \\right )}}\\right|^{2} + \\left|{\\cos{\\left (\\theta \\right )}}\\right|^{2}} \\sqrt{\\left|{p_{T x}}\\right|^{2} + \\left|{p_{T y}}\\right|^{2} + \\left|{p_{T z}}\\right|^{2}}} + \\frac{p_{T z}^{2} \\left(- \\frac{p_{T x} \\cos{\\left (\\theta \\right )}}{\\sqrt{\\left|{\\sin{\\left (\\theta \\right )}}\\right|^{2} + \\left|{\\cos{\\left (\\theta \\right )}}\\right|^{2}} \\sqrt{\\left|{p_{T x}}\\right|^{2} + \\left|{p_{T y}}\\right|^{2} + \\left|{p_{T z}}\\right|^{2}}} - \\frac{p_{T y} \\sin{\\left (\\theta \\right )}}{\\sqrt{\\left|{\\sin{\\left (\\theta \\right )}}\\right|^{2} + \\left|{\\cos{\\left (\\theta \\right )}}\\right|^{2}} \\sqrt{\\left|{p_{T x}}\\right|^{2} + \\left|{p_{T y}}\\right|^{2} + \\left|{p_{T z}}\\right|^{2}}} + 1\\right) \\sin^{2}{\\left (\\theta \\right )}}{\\left(\\left|{\\sin{\\left (\\theta \\right )}}\\right|^{2} + \\left|{\\cos{\\left (\\theta \\right )}}\\right|^{2}\\right) \\left(\\left|{p_{T x}}\\right|^{2} + \\left|{p_{T y}}\\right|^{2} + \\left|{p_{T z}}\\right|^{2}\\right) \\left(\\left|{\\frac{p_{T x} \\sin{\\left (\\theta \\right )}}{\\sqrt{\\left|{\\sin{\\left (\\theta \\right )}}\\right|^{2} + \\left|{\\cos{\\left (\\theta \\right )}}\\right|^{2}} \\sqrt{\\left|{p_{T x}}\\right|^{2} + \\left|{p_{T y}}\\right|^{2} + \\left|{p_{T z}}\\right|^{2}}} - \\frac{p_{T y} \\cos{\\left (\\theta \\right )}}{\\sqrt{\\left|{\\sin{\\left (\\theta \\right )}}\\right|^{2} + \\left|{\\cos{\\left (\\theta \\right )}}\\right|^{2}} \\sqrt{\\left|{p_{T x}}\\right|^{2} + \\left|{p_{T y}}\\right|^{2} + \\left|{p_{T z}}\\right|^{2}}}}\\right|^{2} + \\frac{\\left|{p_{T z} \\sin{\\left (\\theta \\right )}}\\right|^{2}}{\\left(\\left|{\\sin{\\left (\\theta \\right )}}\\right|^{2} + \\left|{\\cos{\\left (\\theta \\right )}}\\right|^{2}\\right) \\left(\\left|{p_{T x}}\\right|^{2} + \\left|{p_{T y}}\\right|^{2} + \\left|{p_{T z}}\\right|^{2}\\right)} + \\frac{\\left|{p_{T z} \\cos{\\left (\\theta \\right )}}\\right|^{2}}{\\left(\\left|{\\sin{\\left (\\theta \\right )}}\\right|^{2} + \\left|{\\cos{\\left (\\theta \\right )}}\\right|^{2}\\right) \\left(\\left|{p_{T x}}\\right|^{2} + \\left|{p_{T y}}\\right|^{2} + \\left|{p_{T z}}\\right|^{2}\\right)}\\right)}\\right)^{2} + sy \\left(- \\frac{p_{T z}^{2} \\left(- \\frac{p_{T x} \\cos{\\left (\\theta \\right )}}{\\sqrt{\\left|{\\sin{\\left (\\theta \\right )}}\\right|^{2} + \\left|{\\cos{\\left (\\theta \\right )}}\\right|^{2}} \\sqrt{\\left|{p_{T x}}\\right|^{2} + \\left|{p_{T y}}\\right|^{2} + \\left|{p_{T z}}\\right|^{2}}} - \\frac{p_{T y} \\sin{\\left (\\theta \\right )}}{\\sqrt{\\left|{\\sin{\\left (\\theta \\right )}}\\right|^{2} + \\left|{\\cos{\\left (\\theta \\right )}}\\right|^{2}} \\sqrt{\\left|{p_{T x}}\\right|^{2} + \\left|{p_{T y}}\\right|^{2} + \\left|{p_{T z}}\\right|^{2}}} + 1\\right) \\sin{\\left (\\theta \\right )} \\cos{\\left (\\theta \\right )}}{\\left(\\left|{\\sin{\\left (\\theta \\right )}}\\right|^{2} + \\left|{\\cos{\\left (\\theta \\right )}}\\right|^{2}\\right) \\left(\\left|{p_{T x}}\\right|^{2} + \\left|{p_{T y}}\\right|^{2} + \\left|{p_{T z}}\\right|^{2}\\right) \\left(\\left|{\\frac{p_{T x} \\sin{\\left (\\theta \\right )}}{\\sqrt{\\left|{\\sin{\\left (\\theta \\right )}}\\right|^{2} + \\left|{\\cos{\\left (\\theta \\right )}}\\right|^{2}} \\sqrt{\\left|{p_{T x}}\\right|^{2} + \\left|{p_{T y}}\\right|^{2} + \\left|{p_{T z}}\\right|^{2}}} - \\frac{p_{T y} \\cos{\\left (\\theta \\right )}}{\\sqrt{\\left|{\\sin{\\left (\\theta \\right )}}\\right|^{2} + \\left|{\\cos{\\left (\\theta \\right )}}\\right|^{2}} \\sqrt{\\left|{p_{T x}}\\right|^{2} + \\left|{p_{T y}}\\right|^{2} + \\left|{p_{T z}}\\right|^{2}}}}\\right|^{2} + \\frac{\\left|{p_{T z} \\sin{\\left (\\theta \\right )}}\\right|^{2}}{\\left(\\left|{\\sin{\\left (\\theta \\right )}}\\right|^{2} + \\left|{\\cos{\\left (\\theta \\right )}}\\right|^{2}\\right) \\left(\\left|{p_{T x}}\\right|^{2} + \\left|{p_{T y}}\\right|^{2} + \\left|{p_{T z}}\\right|^{2}\\right)} + \\frac{\\left|{p_{T z} \\cos{\\left (\\theta \\right )}}\\right|^{2}}{\\left(\\left|{\\sin{\\left (\\theta \\right )}}\\right|^{2} + \\left|{\\cos{\\left (\\theta \\right )}}\\right|^{2}\\right) \\left(\\left|{p_{T x}}\\right|^{2} + \\left|{p_{T y}}\\right|^{2} + \\left|{p_{T z}}\\right|^{2}\\right)}\\right)} - \\frac{\\left(- \\frac{p_{T x} \\sin{\\left (\\theta \\right )}}{\\sqrt{\\left|{\\sin{\\left (\\theta \\right )}}\\right|^{2} + \\left|{\\cos{\\left (\\theta \\right )}}\\right|^{2}} \\sqrt{\\left|{p_{T x}}\\right|^{2} + \\left|{p_{T y}}\\right|^{2} + \\left|{p_{T z}}\\right|^{2}}} + \\frac{p_{T y} \\cos{\\left (\\theta \\right )}}{\\sqrt{\\left|{\\sin{\\left (\\theta \\right )}}\\right|^{2} + \\left|{\\cos{\\left (\\theta \\right )}}\\right|^{2}} \\sqrt{\\left|{p_{T x}}\\right|^{2} + \\left|{p_{T y}}\\right|^{2} + \\left|{p_{T z}}\\right|^{2}}}\\right) \\sqrt{- \\left(\\frac{p_{T x} \\cos{\\left (\\theta \\right )}}{\\sqrt{\\left|{\\sin{\\left (\\theta \\right )}}\\right|^{2} + \\left|{\\cos{\\left (\\theta \\right )}}\\right|^{2}} \\sqrt{\\left|{p_{T x}}\\right|^{2} + \\left|{p_{T y}}\\right|^{2} + \\left|{p_{T z}}\\right|^{2}}} + \\frac{p_{T y} \\sin{\\left (\\theta \\right )}}{\\sqrt{\\left|{\\sin{\\left (\\theta \\right )}}\\right|^{2} + \\left|{\\cos{\\left (\\theta \\right )}}\\right|^{2}} \\sqrt{\\left|{p_{T x}}\\right|^{2} + \\left|{p_{T y}}\\right|^{2} + \\left|{p_{T z}}\\right|^{2}}}\\right)^{2} + 1}}{\\sqrt{\\left|{\\frac{p_{T x} \\sin{\\left (\\theta \\right )}}{\\sqrt{\\left|{\\sin{\\left (\\theta \\right )}}\\right|^{2} + \\left|{\\cos{\\left (\\theta \\right )}}\\right|^{2}} \\sqrt{\\left|{p_{T x}}\\right|^{2} + \\left|{p_{T y}}\\right|^{2} + \\left|{p_{T z}}\\right|^{2}}} - \\frac{p_{T y} \\cos{\\left (\\theta \\right )}}{\\sqrt{\\left|{\\sin{\\left (\\theta \\right )}}\\right|^{2} + \\left|{\\cos{\\left (\\theta \\right )}}\\right|^{2}} \\sqrt{\\left|{p_{T x}}\\right|^{2} + \\left|{p_{T y}}\\right|^{2} + \\left|{p_{T z}}\\right|^{2}}}}\\right|^{2} + \\frac{\\left|{p_{T z} \\sin{\\left (\\theta \\right )}}\\right|^{2}}{\\left(\\left|{\\sin{\\left (\\theta \\right )}}\\right|^{2} + \\left|{\\cos{\\left (\\theta \\right )}}\\right|^{2}\\right) \\left(\\left|{p_{T x}}\\right|^{2} + \\left|{p_{T y}}\\right|^{2} + \\left|{p_{T z}}\\right|^{2}\\right)} + \\frac{\\left|{p_{T z} \\cos{\\left (\\theta \\right )}}\\right|^{2}}{\\left(\\left|{\\sin{\\left (\\theta \\right )}}\\right|^{2} + \\left|{\\cos{\\left (\\theta \\right )}}\\right|^{2}\\right) \\left(\\left|{p_{T x}}\\right|^{2} + \\left|{p_{T y}}\\right|^{2} + \\left|{p_{T z}}\\right|^{2}\\right)}}}\\right)^{2} + sz \\left(\\frac{p_{T z} \\left(- \\frac{p_{T x} \\sin{\\left (\\theta \\right )}}{\\sqrt{\\left|{\\sin{\\left (\\theta \\right )}}\\right|^{2} + \\left|{\\cos{\\left (\\theta \\right )}}\\right|^{2}} \\sqrt{\\left|{p_{T x}}\\right|^{2} + \\left|{p_{T y}}\\right|^{2} + \\left|{p_{T z}}\\right|^{2}}} + \\frac{p_{T y} \\cos{\\left (\\theta \\right )}}{\\sqrt{\\left|{\\sin{\\left (\\theta \\right )}}\\right|^{2} + \\left|{\\cos{\\left (\\theta \\right )}}\\right|^{2}} \\sqrt{\\left|{p_{T x}}\\right|^{2} + \\left|{p_{T y}}\\right|^{2} + \\left|{p_{T z}}\\right|^{2}}}\\right) \\left(- \\frac{p_{T x} \\cos{\\left (\\theta \\right )}}{\\sqrt{\\left|{\\sin{\\left (\\theta \\right )}}\\right|^{2} + \\left|{\\cos{\\left (\\theta \\right )}}\\right|^{2}} \\sqrt{\\left|{p_{T x}}\\right|^{2} + \\left|{p_{T y}}\\right|^{2} + \\left|{p_{T z}}\\right|^{2}}} - \\frac{p_{T y} \\sin{\\left (\\theta \\right )}}{\\sqrt{\\left|{\\sin{\\left (\\theta \\right )}}\\right|^{2} + \\left|{\\cos{\\left (\\theta \\right )}}\\right|^{2}} \\sqrt{\\left|{p_{T x}}\\right|^{2} + \\left|{p_{T y}}\\right|^{2} + \\left|{p_{T z}}\\right|^{2}}} + 1\\right) \\sin{\\left (\\theta \\right )}}{\\sqrt{\\left|{\\sin{\\left (\\theta \\right )}}\\right|^{2} + \\left|{\\cos{\\left (\\theta \\right )}}\\right|^{2}} \\sqrt{\\left|{p_{T x}}\\right|^{2} + \\left|{p_{T y}}\\right|^{2} + \\left|{p_{T z}}\\right|^{2}} \\left(\\left|{\\frac{p_{T x} \\sin{\\left (\\theta \\right )}}{\\sqrt{\\left|{\\sin{\\left (\\theta \\right )}}\\right|^{2} + \\left|{\\cos{\\left (\\theta \\right )}}\\right|^{2}} \\sqrt{\\left|{p_{T x}}\\right|^{2} + \\left|{p_{T y}}\\right|^{2} + \\left|{p_{T z}}\\right|^{2}}} - \\frac{p_{T y} \\cos{\\left (\\theta \\right )}}{\\sqrt{\\left|{\\sin{\\left (\\theta \\right )}}\\right|^{2} + \\left|{\\cos{\\left (\\theta \\right )}}\\right|^{2}} \\sqrt{\\left|{p_{T x}}\\right|^{2} + \\left|{p_{T y}}\\right|^{2} + \\left|{p_{T z}}\\right|^{2}}}}\\right|^{2} + \\frac{\\left|{p_{T z} \\sin{\\left (\\theta \\right )}}\\right|^{2}}{\\left(\\left|{\\sin{\\left (\\theta \\right )}}\\right|^{2} + \\left|{\\cos{\\left (\\theta \\right )}}\\right|^{2}\\right) \\left(\\left|{p_{T x}}\\right|^{2} + \\left|{p_{T y}}\\right|^{2} + \\left|{p_{T z}}\\right|^{2}\\right)} + \\frac{\\left|{p_{T z} \\cos{\\left (\\theta \\right )}}\\right|^{2}}{\\left(\\left|{\\sin{\\left (\\theta \\right )}}\\right|^{2} + \\left|{\\cos{\\left (\\theta \\right )}}\\right|^{2}\\right) \\left(\\left|{p_{T x}}\\right|^{2} + \\left|{p_{T y}}\\right|^{2} + \\left|{p_{T z}}\\right|^{2}\\right)}\\right)} - \\frac{p_{T z} \\sqrt{- \\left(\\frac{p_{T x} \\cos{\\left (\\theta \\right )}}{\\sqrt{\\left|{\\sin{\\left (\\theta \\right )}}\\right|^{2} + \\left|{\\cos{\\left (\\theta \\right )}}\\right|^{2}} \\sqrt{\\left|{p_{T x}}\\right|^{2} + \\left|{p_{T y}}\\right|^{2} + \\left|{p_{T z}}\\right|^{2}}} + \\frac{p_{T y} \\sin{\\left (\\theta \\right )}}{\\sqrt{\\left|{\\sin{\\left (\\theta \\right )}}\\right|^{2} + \\left|{\\cos{\\left (\\theta \\right )}}\\right|^{2}} \\sqrt{\\left|{p_{T x}}\\right|^{2} + \\left|{p_{T y}}\\right|^{2} + \\left|{p_{T z}}\\right|^{2}}}\\right)^{2} + 1} \\cos{\\left (\\theta \\right )}}{\\sqrt{\\left|{\\sin{\\left (\\theta \\right )}}\\right|^{2} + \\left|{\\cos{\\left (\\theta \\right )}}\\right|^{2}} \\sqrt{\\left|{p_{T x}}\\right|^{2} + \\left|{p_{T y}}\\right|^{2} + \\left|{p_{T z}}\\right|^{2}} \\sqrt{\\left|{\\frac{p_{T x} \\sin{\\left (\\theta \\right )}}{\\sqrt{\\left|{\\sin{\\left (\\theta \\right )}}\\right|^{2} + \\left|{\\cos{\\left (\\theta \\right )}}\\right|^{2}} \\sqrt{\\left|{p_{T x}}\\right|^{2} + \\left|{p_{T y}}\\right|^{2} + \\left|{p_{T z}}\\right|^{2}}} - \\frac{p_{T y} \\cos{\\left (\\theta \\right )}}{\\sqrt{\\left|{\\sin{\\left (\\theta \\right )}}\\right|^{2} + \\left|{\\cos{\\left (\\theta \\right )}}\\right|^{2}} \\sqrt{\\left|{p_{T x}}\\right|^{2} + \\left|{p_{T y}}\\right|^{2} + \\left|{p_{T z}}\\right|^{2}}}}\\right|^{2} + \\frac{\\left|{p_{T z} \\sin{\\left (\\theta \\right )}}\\right|^{2}}{\\left(\\left|{\\sin{\\left (\\theta \\right )}}\\right|^{2} + \\left|{\\cos{\\left (\\theta \\right )}}\\right|^{2}\\right) \\left(\\left|{p_{T x}}\\right|^{2} + \\left|{p_{T y}}\\right|^{2} + \\left|{p_{T z}}\\right|^{2}\\right)} + \\frac{\\left|{p_{T z} \\cos{\\left (\\theta \\right )}}\\right|^{2}}{\\left(\\left|{\\sin{\\left (\\theta \\right )}}\\right|^{2} + \\left|{\\cos{\\left (\\theta \\right )}}\\right|^{2}\\right) \\left(\\left|{p_{T x}}\\right|^{2} + \\left|{p_{T y}}\\right|^{2} + \\left|{p_{T z}}\\right|^{2}\\right)}}}\\right)^{2}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(2, 1)\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkEAAAGiCAYAAADgPBIcAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAA9hAAAPYQGoP6dpAABCCUlEQVR4nO3deXRU9f3/8dcQkiEBEiCBLCwhIFQxLhWUxWqIGBYBRZSK9KegCC5EpYALxSVQQFHcvlBcWhpcimBbsdKqLAooBTVGkE0RkEUggRIgAQLJkNzfH9cMxAQIOJPPZO7zcU7OnblzM3m/uQO88rmfe6/LsixLAAAADlPLdAEAAAAmEIIAAIAjEYIAAIAjEYIAAIAjEYIAAIAjEYIAAIAjEYIAAIAjEYIAAIAjEYIAAIAjEYIAAIAj+TUEffrpp+rbt68SEhLkcrn03nvvlXvdsixlZGQoISFB4eHh6tq1q9avX19um6KiIt1///2KiYlR3bp1df3112vnzp3+LBsAADiAX0PQkSNHdMkll2j69OmVvv7MM8/o+eef1/Tp05WVlaW4uDilpaXp0KFD3m1GjhypefPmac6cOVq+fLkOHz6sPn36qKSkxJ+lAwCAIOeqrhuoulwuzZs3T/369ZNkjwIlJCRo5MiReuSRRyTZoz6xsbGaMmWK7r77buXn56tx48Z68803dcstt0iSdu/erebNm+uDDz5Qjx49qqN0AAAQhGqb+sFbt25Vbm6uunfv7l3ndruVkpKiFStW6O6771Z2drY8Hk+5bRISEpScnKwVK1acMgQVFRWpqKjI+7y0tFT79+9XdHS0XC6X/5oCAAA+Y1mWDh06pISEBNWq5fuDV8ZCUG5uriQpNja23PrY2Fht377du01YWJgaNmxYYZuy76/MU089pfHjx/u4YgAAYMKPP/6oZs2a+fx9jYWgMj8fmbEs64yjNWfaZuzYsRo1apT3eX5+vlq0aKHvv/9ejRo1+mUF1yAej0dLlixRamqqQkNDTZdTbeibvp2AvunbCfbv36+2bduqfv36fnl/YyEoLi5Okj3aEx8f712/d+9e7+hQXFyciouLdeDAgXKjQXv37lWXLl1O+d5ut1tut7vC+kaNGik6OtpXLQQ8j8ejiIgIRUdHO+ovDX3TtxPQN307ib+mshi7TlBSUpLi4uK0aNEi77ri4mItW7bMG3Dat2+v0NDQctvk5ORo3bp1pw1BAAAAZ+LXkaDDhw9r8+bN3udbt27V6tWr1ahRI7Vo0UIjR47U5MmT1aZNG7Vp00aTJ09WRESEBg0aJEmKiorS0KFDNXr0aEVHR6tRo0YaM2aMLrroIl177bX+LB0AAAQ5v4agr776Sqmpqd7nZfN0Bg8erFmzZunhhx/W0aNHdd999+nAgQPq2LGjFi5cWO7Y3wsvvKDatWvrt7/9rY4ePapu3bpp1qxZCgkJ8WfpAAAgyPk1BHXt2lWnuwyRy+VSRkaGMjIyTrlNnTp1NG3aNE2bNs0PFQIAAKfi3mEAAMCRCEEAAMCRCEEAAMCRCEEAAMCRCEEAAMCRCEEAAMCRCEEAAMCRCEEAAMCRCEEAAMCRCEEAAMCRCEEAAMCRCEEAAMCRCEEAAMCRCEEAAMCRCEEAAMCRCEEAAMCRCEEAAMCRCEEAAMCRCEEAAMCRCEEAAMCRCEEAAMCRCEEAAMCRCEEAAMCRCEEAAMCRCEEAAMCRCEEAAMCRCEEAAMCRCEEAAMCRCEEAAMCRCEEAAMCRCEEAAMCRCEEAAMCRCEEAAMCRCEEAAMCRCEEAAMCRjIegli1byuVyVfgaMWKEJGnIkCEVXuvUqZPhqgEAQE1X23QBWVlZKikp8T5ft26d0tLSNGDAAO+6nj17KjMz0/s8LCysWmsEAADBx3gIaty4cbnnTz/9tFq3bq2UlBTvOrfbrbi4uOouDQAABDHjIehkxcXFeuuttzRq1Ci5XC7v+qVLl6pJkyZq0KCBUlJSNGnSJDVp0uSU71NUVKSioiLv84KCAkmSx+ORx+PxXwMBpqxXJ/Us0Td9OwN907cT+Ltfl2VZll9/wll45513NGjQIO3YsUMJCQmSpLlz56pevXpKTEzU1q1b9fjjj+v48ePKzs6W2+2u9H0yMjI0fvz4Cutnz56tiIgIv/YAAAB8o7CwUIMGDVJ+fr4iIyN9/v4BFYJ69OihsLAwzZ8//5Tb5OTkKDExUXPmzFH//v0r3aaykaDmzZsrJydH0dHRPq87UHk8Hi1atEhpaWkKDQ01XU61oW/6dgL6pm8nyMvLU3x8vN9CUMAcDtu+fbsWL16sd99997TbxcfHKzExUZs2bTrlNm63u9JRotDQUEd9eMrQt7PQt7PQt7M4rW9/92r8FPkymZmZatKkiXr37n3a7fLy8vTjjz8qPj6+mioDAADBKCBCUGlpqTIzMzV48GDVrn1icOrw4cMaM2aMVq5cqW3btmnp0qXq27evYmJidOONNxqsGAAA1HQBcThs8eLF2rFjh+68885y60NCQrR27Vq98cYbOnjwoOLj45Wamqq5c+eqfv36hqoFAADBICBCUPfu3VXZ/Ozw8HAtWLDAQEUAACDYBcThMAAAgOpGCAIAAI5ECAIAAI5ECAIAAI5ECAIAAI5ECAIAAI5ECAIAAI5ECAIAAI5ECAIAAI5ECAIAAI5ECAIAAI4UEPcOA4BzVloq94EDcs2fL9WuLV14odSwoRQVJYWEmK4OQAAjBAGoOY4elRYulP71L+nbb6Vdu1Q7J0c9jx+vuK3LJTVoILVqJXXvLv3xj4QiAOUQggDUHMuWSf36lVvlkmS5XHJZlr2iXj3p8GHJsqQDB6TsbKmoSJo8+cQ3ffihdOWVUmRktZUOIPAQggAEpl27pPHjpSZNpIkT7XXXXCNdfLG9/M1vpGbN5GnSRB+uWqVeffsqNDTU3q642A5AeXnS11+XHwE6elTq318KD5dGj5buv58wBDgUIQhAYDl4UJoyRXrpJTuwtG5tH8pyuaSwMOmbb8pv7/HIWrOm/LqwMCk21v5q1678azk5UmKitHGj9Nhj0nPPEYYAh+LsMACB4dgxO5C0aiU9/bQdgK68Unr9dTsA+UqrVtL69dLs2dL559sjRo89JrVsKU2dKpWU+O5nAQhohCAA5q1ZYx/mGjPGDiXt2knvvy999pkdhHwtJES69VZp3bryYeiRR6TVq33/8wAEJA6HATDv0CFp61YpIcGe/3P77dVzJldZGPrtb6XMTGnfPql9e///XAABgRAEwAzLOnGY68orpb//XbrqKik6uvprCQmR7rqr/Lrvv5f+/Gc7lLnd1V8TAL/jcBiA6vfDD1JKij03p0y/fmYCUGVKSuwRoqlTpc6d7cnUAIIOIQhA9crOljp0sOf73HOPPSIUaEJC7NPzY2KkVavswLZzp+mqAPgYIQhA9Vm1SkpLsychX3GFPSnZl2d++VKfPtIXX9in02/aZAeh7dtNVwXAhwhBAKrHN99I115rB6DOnaXFi6XmzU1XdXqtWtlXqW7Vyj6Ed/XV9hJAUCAEAfC/devsALR/vz0C9NFHUv36pquqmsRE6dNPpbZtpR077IsqAggKhCAA/jdu3InTzxcsqHlXZm7aVFq6VLrpJvtUegBBgRAEwP/efFMaPty+A3yDBqarOTfx8dI//mHfywxAUCAEAfC/yEjp1VelRo1MV+I7f/mLNGRIYJ7dBqBKuFgiAP8YNcq+AvTo0YF7Bti52rxZuvde6fhx6cILpYceMl0RgHPASBAA33v7bemFF+xw8OWXpqvxvfPOs+9yL0mPPmof5gNQ4xCCAPjW+vUnbkHxhz9IHTuarcdf7r1XGjpUKi2VBg6UtmwxXRGAs0QIAuA7hw7ZZ1AVFkrdukkTJpiuyH9cLulPf5I6dbKvfXTjjVJRkemqAJwFQhAA3xk+XNq4UWrWzD4kVh13gjfJ7Zb++U/7jLG1a6VJk0xXBOAsEIIA+Mb8+dKcOXbweecdqXFj0xVVj4QEe0QoJIQzxYAahrPDAPjGtm1SWJg0cqR9Wwwnuekm6dtvpTZtTFcC4CwYHwnKyMiQy+Uq9xUXF+d93bIsZWRkKCEhQeHh4eratavWr19vsGIAlbr/fmnNGumJJ0xXUv1cLgIQUAMZD0GSdOGFFyonJ8f7tXbtWu9rzzzzjJ5//nlNnz5dWVlZiouLU1pamg4dOmSwYgCV+tWvpLp1TVdh1vr19qRwzhYDAl5AhKDatWsrLi7O+9X4p7kElmXpxRdf1Lhx49S/f38lJyfr9ddfV2FhoWbPnm24agAqKbFPE8/KMl1J4BgzRvrkE3uSOHOEgIAWEHOCNm3apISEBLndbnXs2FGTJ09Wq1attHXrVuXm5qp79+7ebd1ut1JSUrRixQrdfffdlb5fUVGRik46VbWgoECS5PF45PF4/NtMACnr1Uk9S/RdnX3Xevllhfz1r7LmzdPxLVukevWq7WeXCbj9/cILqr1smVyffKLjmZmybrvNLz8m4PquJvTtzL79xWVZZn9V+fDDD1VYWKi2bdtqz549mjhxor777jutX79eGzdu1JVXXqldu3YpISHB+z3Dhw/X9u3btWDBgkrfMyMjQ+PHj6+wfvbs2YqIiPBbL4CTuPfvV7cRIxR69Ki+GT5c2667znRJAeO8d9/VhW+8ocLGjfXxjBkqDQ01XRJQIxUWFmrQoEHKz89XZGSkz9/feAj6uSNHjqh169Z6+OGH1alTJ1155ZXavXu34uPjvdsMGzZMP/74oz766KNK36OykaDmzZsrJydH0dHRfu8hUHg8Hi1atEhpaWkKddA/wvRdPX3XevBBhbz8sko7dFDJZ58ZuyZQQO7vo0dVu107uXbtUslzz6n0/vt9/iMCsu9qQN/O6jsvL0/x8fF+C0EBcTjsZHXr1tVFF12kTZs2qV+/fpKk3NzcciFo7969io2NPeV7uN1uud3uCutDQ0Md9eEpQ9/OUi1979ghzZwpSar1zDOqVaeOf39eFQTU/g4NlZ58Uho+XCFPP62QYcOk+vX99KMCqO9qRN/O4O9eA2Ji9MmKior07bffKj4+XklJSYqLi9OiRYu8rxcXF2vZsmXq0qWLwSoBh5s8WSoullJT7S9UdMcd9mnz//ufNzACCCzGR4LGjBmjvn37qkWLFtq7d68mTpyogoICDR48WC6XSyNHjtTkyZPVpk0btWnTRpMnT1ZERIQGDRpkunTAmbZtO/GfeiVz7/CT2rWl55+Xdu+2AxGAgGM8BO3cuVO33nqr9u3bp8aNG6tTp076/PPPlZiYKEl6+OGHdfToUd133306cOCAOnbsqIULF6q+n4aWAZxB06bSK69IK1ZIV11luprA1qeP6QoAnIbxEDRnzpzTvu5yuZSRkaGMjIzqKQjA6YWG2tcGGjrUdCU1S3GxdPSoFBVluhIAPwm4OUEAAlhpqekKaqaFC6Xzz5fGjTNdCYCTEIIAVM3GjfZE35kzuRLy2QoJkbZulWbNkg4eNF0NgJ8QggBUzUsvST/8IP3rX/YNQ1F111wjJSdLR45If/2r6WoA/IQQBODMCgqkN9+0H48cabSUGsnlkh54wH48fbp9zzUAxhGCAJzZm29Khw9LF1zAdYHO1e9+JzVsaB8W+89/TFcDQIQgAGdiWdKf/mQ/vu8+DoWdq4gIadgw+/FLL5mtBYAkQhCAM1m2TPr2W6luXclPd0R3jBEjpFq1pE8+sSeaAzDK+HWCAAS4slGg227jGje/VIsW0rPPSldcIbVta7oawPEIQQBO77777OsD3Xef6UqCw6hRpisA8BNCEIDT4yapAIIUc4IAoLrt3Cmlp0vXX2+6EsDRCEEAKvff/0pjx0pr15quJPiEhEgvvyzNn29fgBKAEYQgAJXLzJSeftq+uB98Kz7evoq0JL39ttlaAAcjBAGoqLhYevdd+/Ett5itJVgNGmQv//Y37sUGGEIIAlDRxx9LBw5IsbFSSorpaoJT//6S221fg2nNGtPVAI5ECAJQ0dy59vLmm+35K/C9qCipd2/78ezZZmsBHIoQBKC8oiLpvffsxxwK86+yQ2Jvv21fiwlAteI6QQDKW7BAys+XmjaVrrzSdDXBrXdvKTlZuu466ehR+9YkAKoNIQhAefv2SdHR0oAB9n2u4D916nAJAsAgQhCA8u68075PWGGh6UoAwK/4NQ9ARaGh3Cy1OhUVSYsXS7t2ma4EcBRCEIAT9u7lmjUm3HSTlJYmvfOO6UoARyEEATjh6qulZs2k7GzTlThL2dWjP/rIbB2AwxCCANh27JA2bpRyc6XWrU1X4yy9etnLZcuYiwVUI0IQANuiRfayY0epQQOjpTjO+edLLVrYc4OWLjVdDeAYhCAAtoUL7WVamtk6nMjlOjEaxCExoNoQggBIJSX22UmS1L272VqcqmdPe/nhh2brAByEEARAWrVK2r9fioyUrrjCdDXO1K2bVLu2tHmztGWL6WoAR+BiiQBOHApLTbWvEYTqV7++NGuWdNFFUqtWpqsBHIEQBEC69lp7JKhjR9OVONvvfme6AsBRCEEA7ENgHAYD4DDMCQKAQPLXv0pDh0p79piuBAh6hCDA6b780r5G0MGDpiuBJL3wgh2EVqwwXQkQ9AhBgNO9+KJ9Wvz06aYrgSR17mwvV640WwfgAIQgwOm++MJeMik6MBCCgGpDCAKcbN8+6Ycf7McdOpitBbayEPTVV1JxsdlagCBnPAQ99dRTuvzyy1W/fn01adJE/fr108aNG8ttM2TIELlcrnJfnTp1MlQxEESysuxl27ZSw4Zma4GtbF8cOyZ9843paoCgZjwELVu2TCNGjNDnn3+uRYsW6fjx4+revbuOHDlSbruePXsqJyfH+/XBBx8YqhgIIhwKCzy1akllv+R9/rnZWoAgZ/w6QR/97GaBmZmZatKkibKzs3X11Vd717vdbsXFxVXpPYuKilRUVOR9XlBQIEnyeDzyeDw+qLpmKOvVST1L9H02fYd88YVqSSpp316lNfTPKxj3d60rrlDIhx+qZMuWU+6XYOy7KujbmX37i8uyLMuvP+Esbd68WW3atNHatWuVnJwsyT4c9t577yksLEwNGjRQSkqKJk2apCZNmlT6HhkZGRo/fnyF9bNnz1ZERIRf6wdqDMtSz9tvl/vQIS179lkdbNPGdEX4SWhBgeRyyVO/vulSAKMKCws1aNAg5efnKzIy0ufvH1AhyLIs3XDDDTpw4IA+++wz7/q5c+eqXr16SkxM1NatW/X444/r+PHjys7OltvtrvA+lY0ENW/eXDk5OYqOjq6WXgKBx+PRokWLlJaWplAH3Q+KvqvYt2VJ33yjWllZKr39dqmSv0s1Afubvp3AqX3n5eUpPj7ebyHI+OGwk6Wnp2vNmjVavnx5ufW33HKL93FycrI6dOigxMRE/ec//1H//v0rvI/b7a40HIWGhjrqw1OGvp3lrPq+/HLp8ssV4t+SqgX721no2xn83avxidFl7r//fr3//vtasmSJmjVrdtpt4+PjlZiYqE2bNlVTdQBQzf7v/6Ru3aT33jNdCRC0jIcgy7KUnp6ud999V5988omSkpLO+D15eXn68ccfFR8fXw0VAkHqjTekGTOkbdtMV4LKbNggffKJfVsTAH5hPASNGDFCb731lmbPnq369esrNzdXubm5Onr0qCTp8OHDGjNmjFauXKlt27Zp6dKl6tu3r2JiYnTjjTcarh6owV56SRoxQlq92nQlqMyFF9rL9evN1gEEMeNzgl5++WVJUteuXcutz8zM1JAhQxQSEqK1a9fqjTfe0MGDBxUfH6/U1FTNnTtX9TlzAjg3paXSd9/Zjy+4wGwtqFy7dvZywwazdQBBzHgIOtPJaeHh4VqwYEE1VQM4xK5dUmGhVLu21KqV6WpQmbKRoC1bpKNHpfBws/UAQcj44TAABnz7rb087zzJQWea1CixsVKjRvalDMpG7QD4FCEIcCIOhQU+l+vEaBCHxAC/IAQBTlQWgs4/32wdOL127aQmTexDlwB8zvicIAAGlB0OIwQFtmnTpFdeMV0FELQIQYAT/eMf9mjQeeeZrgSnw3wtwK8IQYATRUdLV15pugoAMIo5QQAQqAoL7VtntG0rnXRTaAC+QQgCnGb9eunhh6W//c10JTiT8HBp5Upp0yZp507T1QBBhxAEOM3XX0vPPivNmmW6EpyJyyW1aGE/3rHDbC1AECIEAU5TNqLQtKnZOlA1zZvbyx9/NFsHEIQIQYDT7NplLwlBNQMjQYDfEIIApyEE1SyMBAF+QwgCnKYsBDVrZrYOVE1ZCGIkCPA5QhDgNMwJqlkSE+1bZ0RFma4ECDpcLBFwkuPHpT177MeEoJrh2mtP7DMAPkUIApwkJMQ+rLJzpz26AAAORggCnMTlskeAGAUCAOYEAUDA69VLatdO2rLFdCVAUCEEAU7y1VfSI49wy4ya5rvvpG+/lf73P9OVAEGFEAQ4SXa29Mwz0jvvmK4EZ6NhQ3t54IDZOoAgQwgCnKSgwF5yunXN0qCBvTx40GQVQNAhBAFOkp9vLwlBNQsjQYBfEIIAJyEE1UyEIMAvCEGAk3A4rGYiBAF+QQgCnKRsJCgy0mwdODtxcfa1nSIiTFcCBBUulgg4CYfDaqbRo+0vAD5FCAKc5M03pX37pBYtTFcCAMYRggAnadbM/gIAMCcIAALep59KHTtKd9xhuhIgqDASBDjJlClSUZF0771S48amq0FVFRRIX34plZaargQIKoQgwEmmTLFPs/7tbwlBNUloqL08ftxsHUCQ4XAY4CQlJfayFn/1a5TaP/2+6vGYrQMIMowEAU5SdjglJMRsHX7w5db92ld4XE3q19EVSY0UUstluiTfYSQI8AtCEOAkZSNBQRSCFn+7R5J05+tZKiqxg098VB092bedeibHmyzNd8pCECNBgE8xJg44SdlIUJAcDvtoXY5+P3d1hfW5+cd071tf66N1OdVflD+UHQ5jJAjwqRrzL+GMGTOUlJSkOnXqqH379vrss89MlwTUPEE0ElRSamn8/A2yKnmtbN34+RtUUlrZFjVMnTr2/cMaNDBdCRBUakQImjt3rkaOHKlx48Zp1apVuuqqq9SrVy/t2LHDdGlAzRJEI0Ffbt2vnPxjp3zdkpSTf0xfbt1ffUX5y0UXSfv3S998Y7oSIKjUiDlBzz//vIYOHaq77rpLkvTiiy9qwYIFevnll/XUU09V2L6oqEhFRUXe5wU/3Tnb4/HI46Bj6mW9Oqlnib5P17dr+XKptFRWZGSNn1+yN/+I3CGW3LXskZ6yZWXbeTzBd8NYPuf07QT+7tdlWVZAjxUXFxcrIiJCf//733XjjTd61z/44INavXq1li1bVuF7MjIyNH78+ArrZ8+erQjuwgwAQI1QWFioQYMGKT8/X5GRvv9lJuBHgvbt26eSkhLFxsaWWx8bG6vc3NxKv2fs2LEaNWqU93lBQYGaN2+u1NRURUdH+7XeQOLxeLRo0SKlpaUptOzsEgegb2f0XVJqqceLn+rg4aOa0KFUj39VS0WlJ06Ld0mKjayjBSOvrvmny2/apJC775bi41Xyt79Jct7+LkPfzuo7Ly/Pr+8f8CGojMtV/h8xy7IqrCvjdrvldrsrrA8NDXXUh6cMfTvLaft+6SWpsFC6/36pXr3qLczHQiWN7X2hRr6dLUkqKnV5T5Ev+5dhbO8LVccdZqZAXzp8WFq+XGrZUrV+tm/5nDuL0/r2d68BPzsyJiZGISEhFUZ99u7dW2F0CMAZPPKI9Ic/2JNsg0DP5Hi9cMulFdbHRdXRy//vsuC5TtCxnyaAV/LLHYBzF/AjQWFhYWrfvr0WLVpUbk7QokWLdMMNNxisDKiBIiLsG6gWFpquxGeuvSBWH2yV/jr48uC9YnRZCKpTx2wdQJAJ+BAkSaNGjdJtt92mDh06qHPnznrttde0Y8cO3XPPPaZLA2qW8HD7BqpHj5quxOeuSGoUvIcJys52JQQBPlUjQtAtt9yivLw8TZgwQTk5OUpOTtYHH3ygxMRE06UBNUvZ2ZFBNBLkCBwOA/yiRoQgSbrvvvt03333mS4DqNkIQTUTh8MAvwj4idEAfCg83F4G4eGwoGZZdoDlOmeATxGCACdhJKhmuv126cgRad4805UAQaXGHA4D4APPPmv/Z3rBBaYrAQDjCEGAk7Rvb7oCAAgYHA4DgED3f/8nXXedNGeO6UqAoEIIApxkxQpp+nT7FgyoOb7+WvrwQ2nHDtOVAEGFEAQ4yQcf2PcNmzvXdCU4G/v22cuYGLN1AEGGEAQ4SZMm9nLvXrN14OwQggC/IAQBTkIIqpkIQYBfEIIAJyEE1UyEIMAvCEGAkxCCah6PR8rPtx9HR5utBQgyhCDAScpCUF6edPy42VpQNQcPSpGRUq1aUoMGpqsBggohCHCS6GjJ5bLvRZWXZ7oaVEXjxvZI0NGjUkiI6WqAoMIVowEnCQmR5s+XGjZkVKGmCQszXQEQdAhBgNP07m26AgAICBwOA4BANnu21KOH9MorpisBgg4jQYDTrFghffWV1KGD1KWL6WpwJl9/LS1cKF14oelKgKDDSBDgNH//u/Tgg9K8eaYrQVVs3Wovk5LM1gEEIUIQ4DStWtnLLVvM1oGqIQQBfkMIApymdWt7+cMPZutA1ZSFoJYtjZYBBCNCEOA0ZSFoyxb7ekEIXAcP2l8SIQjwA0IQ4DQtW9oXTDx8WPrf/0xXg9PZts1eNm4s1atntBQgGBGCAKdxu6VmzezHzAsKbAcO2Ff5ZhQI8AtCEOBEJx8SQ+BKTbXvIP/ZZ6YrAYIS1wkCnGjKFHt5wQVm60DVuN2mKwCCEiEIcKIrrjBdAQAYx+EwAAhUffpIPXtKGzaYrgQISowEAU5UUCDNmiXt2SNNmmS6GlSmqEj66COppESKjDRdDRCUCEGAE5WW2rfOkKSHH5aioszWg4q+/dYOQA0bSk2bmq4GCEocDgOcqEEDqXlz+/GaNUZLwSmU7ZeLL7av6wTA5whBgFNdeqm9XLXKaBk4hZNDEAC/IAQBTvXrX9vL1auNloFTIAQBfkcIApyKkaDARggC/I6J0YBTlY0ErV8vFRdLYWFm68EJhYXSeefZZ4hdeKHpaoCgxUgQ4FSJifYEaY9H+v5709XgZBER0vLl0v79Ut26pqsBgpaxELRt2zYNHTpUSUlJCg8PV+vWrfXkk0+quLi43HYul6vC1yuvvGKoaiCIuFzSkiX2TTqTk01Xg8pwVhjgV8YOh3333XcqLS3Vq6++qvPOO0/r1q3TsGHDdOTIEU2dOrXctpmZmerZs6f3eRTXNAF8o2xeEAJLSYkUEmK6CiDoGQtBPXv2LBdsWrVqpY0bN+rll1+uEIIaNGiguLi46i4RAMxo3146elR66y3p8stNVwMErYCaGJ2fn69GjRpVWJ+enq677rpLSUlJGjp0qIYPH65atU59JK+oqEhFRUXe5wUFBZIkj8cjj8fj+8IDVFmvTupZou+z6js/X7UmTpTr++9V8t57NfLwS9Dt7wMHVHvNGrksS574eHvOViWCru8qom9n9u0vLsuyLL/+hCrasmWLLrvsMj333HO66667vOsnTpyobt26KTw8XB9//LGeeOIJjR07Vo899tgp3ysjI0Pjx4+vsH727NmKiIjwS/1ATeQ6fly9Bw1SSHGxPp42TYfLriINY2KzstRp0iQdTkjQxzNmmC4HMKqwsFCDBg1Sfn6+Iv1wDz2fh6BTBZCTZWVlqUOHDt7nu3fvVkpKilJSUvSXv/zltN/73HPPacKECcrPzz/lNpWNBDVv3lw5OTmKjo6uYic1n8fj0aJFi5SWlqbQ0FDT5VQb+j67vkPS0lRr2TKVTJ+u0uHD/VihfwTb/q71hz8oZOpUld5xh0peffWU2wVb31VF387qOy8vT/Hx8X4LQT4/HJaenq6BAweedpuWLVt6H+/evVupqanq3LmzXnvttTO+f6dOnVRQUKA9e/YoNja20m3cbrfcbneF9aGhoY768JShb2c56767dpWWLVPI8uUKGTHCb3X5W9Ds7xUrJEm1rr5atarQT9D0fZbo2xn83avPQ1BMTIxiYmKqtO2uXbuUmpqq9u3bKzMz87TzfMqsWrVKderUUYMGDX5hpQAkSSkp9nLZMsmyauS8oKBx9KiUlWU/vuoqs7UADmBsYvTu3bvVtWtXtWjRQlOnTtX//vc/72tlZ4LNnz9fubm56ty5s8LDw7VkyRKNGzdOw4cPr3SkB8A56NTJvlp0To60ebPUpo3pipwrK8ueCB0XJ7VqZboaIOgZC0ELFy7U5s2btXnzZjVr1qzca2XTlEJDQzVjxgyNGjVKpaWlatWqlSZMmKARNXjIHgg44eHSFVfYVyhetowQZFJUlDR0qL1kRA7wO2MhaMiQIRoyZMhpt/n5tYQA+EnXrtLGjfbhGJhzySXSGU4OAeA7AXWdIACGjBsnTZjA6AMAR+EGqgCkOnUIQKbl5EhffSUdP266EsAxCEEATrAs6cgR01U4U2amfYuMW24xXQngGIQgALY5c6SEBOmee0xX4kz//re97N7dbB2AgxCCANgSEqTcXOmDDzgkU9327pU+/9x+3Lu32VoAByEEAbB16SJFR0v799uny6P6fPihfSjy17+WfnbJEAD+QwgCYKtdW+rb13783ntGS3Gc+fPtZZ8+ZusAHIYQBOCEG26wl//6lz0yAf8rLpYWLLAfl4VQANWCEATghO7d7StIb9smrVljuhpn+PRT6fBh+1YZ7dubrgZwFEIQgBMiIqS0NPvxv/5lthanSE2V/vtfafp0qQo3kQbgO1wxGkB5d9xh3z+M+SnVIyTEnpQOoNoRggCU16+f/QUAQY6xVwAwZeJEadgwad0605UAjkQIAlBRaan0ySfS+PGmKwlex49LM2bYd43fvNl0NYAjcTgMQEV79thnipWUSAMHSr/6lemKgs/ixfZNU6OjpeuuM10N4EiMBAGoKD5e6tXLfpyZabaWYPXGG/by1lulsDCztQAORQgCULk777SXr7/OvcR8LT9fmjfPfjx4sNlaAAcjBAGoXO/eUuPG9k1VP/rIdDXB5e9/l44dky64gAskAgYRggBULixMuu02+/Ff/2q2lmDz+uv2cvBgyeUyWwvgYIQgAKd2xx32cv58ae9es7UEi9JS6ZprpFatpN/9znQ1gKMRggCcWnKydMUV9n/Y27ebriY41KplX3pg82apWTPT1QCOxinyAE7v3/+WYmI4bONr/HkCxjESBOD0GjfmP2xf+fe/pX/+077+EgDjCEEAquboUe4s/0uUlkoPPSTdfLP05z+brgaAOBwGoCqOHLHvLJ+TI61aJV16qemKap4PPpC++06KimJCNBAgGAkCcGZ160pdu9qPp0wxWkqNNXWqvbz7bql+fbO1AJBECAJQVY88Yi/feUfassVsLTVNVpa0bJlUu7b0wAOmqwHwE0IQgKq55BL7fmKlpSdGNVA1ZX9egwZJTZuarQWAFyEIQNU9+qi9zMy0b6eBM9u6VfrHP+zHo0ebrQVAOYQgAFV31VVS585SUZH00kumq6kZjh6VLr9cSkuTLr7YdDUATkIIAlB1LteJ0aD8fLO11BTt2kkrV0pvv226EgA/wynyAM5Onz7SV19x9/Oz4XJJ0dGmqwDwM4wEATg7tWoRgKpi8WLp8celQ4dMVwLgFAhBAM7dtm3SxImSZZmuJLCUlEijRtl/Nk8/bboaAKfA4TAA5+bwYenXv5YOHpQuuki64QbTFQWO11+X1q6VGjbkjDAggBkdCWrZsqVcLle5r0fLJl3+ZMeOHerbt6/q1q2rmJgYPfDAAyouLjZUMQCvevWke++1H48ebZ8xBvsWI48/bj9+7DGpUSOz9QA4JeMjQRMmTNCwYcO8z+vVq+d9XFJSot69e6tx48Zavny58vLyNHjwYFmWpWnTppkoF8DJxo61rxm0ZYs0bZo0Zozpisx75hlp924pKUkaMcJ0NQBOw/icoPr16ysuLs77dXIIWrhwoTZs2KC33npLv/71r3Xttdfqueee05///GcVFBQYrBqAJPseWJMn24//+EcuoLhq1Yk/jylTJLfbbD0ATsv4SNCUKVP0xz/+Uc2bN9eAAQP00EMPKSwsTJK0cuVKJScnKyEhwbt9jx49VFRUpOzsbKWmplb6nkVFRSo6aWi+LDB5PB55PB4/dhNYynp1Us8SfVd734MGKWT6dNX6+muVDh2qknnz7FPCq0kg7e+QkSNV6/hxld54o0puuEHyY02B1Hd1om9n9u0vLssyd1rHCy+8oMsuu0wNGzbUl19+qbFjx+qGG27QX/7yF0nS8OHDtW3bNi1cuLDc97ndbs2aNUu33nprpe+bkZGh8ePHV1g/e/ZsRURE+L4RwOHqb9umlIceUojHo69GjdKuq682XZIR7v37deEbb2jdHXeoOCrKdDlAjVdYWKhBgwYpPz9fkZGRPn9/n4egUwWQk2VlZalDhw4V1v/zn//UzTffrH379ik6OlrDhw/X9u3btWDBgnLbhYWF6Y033tDAgQMrff/KRoKaN2+unJwcRTvogmUej0eLFi1SWlqaQkNDTZdTbejbTN+1/vQnKSdHpU88If00mlsdTPdtCn3TtxPk5eUpPj7ebyHI54fD0tPTTxlOyrRs2bLS9Z06dZIkbd68WdHR0YqLi9MXX3xRbpsDBw7I4/EoNjb2lO/vdrvlruRYfGhoqKM+PGXo21mM9T1ypCQppPp/siSDfR85Ii1bJl13XfX/bPE5dxqn9e3vXn0egmJiYhQTE3NO37tq1SpJUnx8vCSpc+fOmjRpknJycrzrFi5cKLfbrfZcsRYIXB6PtGSJ1L276Ur879FHpenTpYcftidDA6gxjE2MXrlypT7//HOlpqYqKipKWVlZ+v3vf6/rr79eLVq0kCR1795d7dq102233aZnn31W+/fv15gxYzRs2DC/DIsB8IHCQiklRcrOtoNQSorpivznH/+wA5AkdetmthYAZ83YKfJut1tz585V165d1a5dOz3xxBMaNmyY3j7pTsshISH6z3/+ozp16ujKK6/Ub3/7W/Xr109Tp041VTaAM4mIkC6+2L6Vxq23Sjt2mK7IP778UrrtNvvx73/vjFEvIMgYGwm67LLL9Pnnn59xuxYtWujf//53NVQEwGdeekn64gtp/XqpVy9p+XL7FhLBYvt26frrpWPHpN69pWefNV0RgHNg/GKJAIJQvXrShx9KCQnShg1Sv352YAgGBQVSnz7Snj32iNfbb0shpqaDA/glCEEA/KN5czsIRUZKn34qDR4slZaaruqXW7hQWrdOio+X/v1v+6rZAGokQhAA/7n4YmnePCk0VFq8WNq61XRFv9zNN0v//Kc0f74d9ADUWMZvmwEgyF1zjTR3rtSundS6telqzs3x49K+fVJcnP28f3+z9QDwCUaCAPjfjTdKv/rViedffy2VlJir52wcPSoNGCBdfbW0d6/pagD4ECEIQPX673+lLl2km26yrykUyA4elHr2lN57zz7Vf80a0xUB8CFCEIDqtWePvfzXv6SuXU88DzQ5OfaFHj/91J7cvWCBdO21pqsC4EOEIADVq39/e5J0o0ZSVpbUubP03Xemqyrv00/tutaskWJj7efBfOVrwKEIQQCq329+I61caU+U3rpVat9emjrVnoBs2ty5duDZvl067zxpxQrpkktMVwXADwhBAMxo29YOQtdcY88Neugh6f33TVdlzwFq3ly66y77/metWpmuCICfcIo8AHMaN7YPjWVm2nNubrzxxGuWJblc/q/hwAHpnXek4cPtnxcVZR8Ga9DA/z8bgFGMBAEwy+WS7rzTPgxVFnry86VOnaTJk/03cXrLFumBB+xRn3vukWbOPPEaAQhwBEIQgMDzzDP2XdrHjbNDyi23SEuX2qNDv4Rl2Tdz7d9fatNGmjZNOnJEuugiqVkzn5QOoOYgBAEIPI89Jr3+uj0a5PHYh6tSU+2rTv/+9/YhrLNVUiKdf7501VX2rTwsy77D/aJF0jff2HOBADgKc4IABJ7wcOn22+2v1aulV16R3nrLPpV+82bp6ae9m9aaOlXtP/xQIW+9ZQed48ft4FRYeOKeZZJ9p/cGDSS3W7rtNjtMtWtnpD0AgYEQBCCwXXqpHYKeeca+wOKuXXaQ+Ynrn/9Us+zsyr83JEQ6fFiqV89+/uab9v2/IiP9XzeAgEcIAlAzREbaIzg/U/rQQ1q/cKHaXXSRQurUkWrXtkeAQkOlCy6wR5XKtG1bjQUDCHSEIAA1mtW/v36oU0fnX3edQkJDTZcDoAZhYjQAAHAkQhAAAHAkQhAAAHAkQhAAAHAkQhAAAHAkQhAAAHAkQhAAAHAkQhAAAHAkQhAAAHAkQhAAAHAkQhAAAHAkQhAAAHAkQhAAAHAkQhAAAHAkQhAAAHAkQhAAAHAkQhAAAHAkYyFo6dKlcrlclX5lZWV5t6vs9VdeecVU2QAAIEjUNvWDu3TpopycnHLrHn/8cS1evFgdOnQotz4zM1M9e/b0Po+KiqqWGgEAQPAyFoLCwsIUFxfnfe7xePT+++8rPT1dLper3LYNGjQoty0AAMAvZSwE/dz777+vffv2aciQIRVeS09P11133aWkpCQNHTpUw4cPV61apz6SV1RUpKKiIu/zgoICSXbQ8ng8Pq89UJX16qSeJfqmb2egb/p2An/367Isy/LrT6ii6667TpL0wQcflFs/ceJEdevWTeHh4fr444/1xBNPaOzYsXrsscdO+V4ZGRkaP358hfWzZ89WRESEbwsHAAB+UVhYqEGDBik/P1+RkZE+f3+fh6BTBZCTZWVllZv3s3PnTiUmJuqdd97RTTfddNrvfe655zRhwgTl5+efcpvKRoKaN2+unJwcRUdHV7GTms/j8WjRokVKS0tTaGio6XKqDX3TtxPQN307QV5enuLj4/0Wgnx+OCw9PV0DBw487TYtW7Ys9zwzM1PR0dG6/vrrz/j+nTp1UkFBgfbs2aPY2NhKt3G73XK73RXWh4aGOurDU4a+nYW+nYW+ncVpffu7V5+HoJiYGMXExFR5e8uylJmZqdtvv71Kza5atUp16tRRgwYNfkGVAADA6YxPjP7kk0+0detWDR06tMJr8+fPV25urjp37qzw8HAtWbJE48aN0/Dhwysd6QEAAKgq4yFo5syZ6tKliy644IIKr4WGhmrGjBkaNWqUSktL1apVK02YMEEjRowwUCkAAAgmxkPQ7NmzT/laz549y10kEQAAwFe4dxgAAHAkQhAAAHAkQhAAAHAkQhAAAHAkQhAAAHAkQhAAAHAkQhAAAHAkQhAAAHAkQhAAAHAkQhAAAHAkQhAAAHAkQhAAAHAkQhAAAHAkQhAAAHAkQhAAAHAkQhAAAHAkQhAAAHAkQhAAAHAkQhAAAHAkQhAAAHAkQhAAAHAkQhAAAHAkQhAAAHAkQhAAAHAkQhAAAHAkQhAAAHAkQhAAAHAkQhAAAHAkQhAAAHAkQhAAAHAkQhAAAHAkQhAAAHAkQhAAAHAkQhAAAHAkQhAAAHAkQhAAAHAkv4agSZMmqUuXLoqIiFCDBg0q3WbHjh3q27ev6tatq5iYGD3wwAMqLi4ut83atWuVkpKi8PBwNW3aVBMmTJBlWf4sHQAABLna/nzz4uJiDRgwQJ07d9bMmTMrvF5SUqLevXurcePGWr58ufLy8jR48GBZlqVp06ZJkgoKCpSWlqbU1FRlZWXp+++/15AhQ1S3bl2NHj3an+UDAIAg5tcQNH78eEnSrFmzKn194cKF2rBhg3788UclJCRIkp577jkNGTJEkyZNUmRkpP72t7/p2LFjmjVrltxut5KTk/X999/r+eef16hRo+RyuSq8b1FRkYqKirzP8/PzJUn79+/3cYeBzePxqLCwUHl5eQoNDTVdTrWhb/p2Avqmbyco+3/bX0d//BqCzmTlypVKTk72BiBJ6tGjh4qKipSdna3U1FStXLlSKSkpcrvd5bYZO3astm3bpqSkpArv+9RTT3kD2Mnatm3rn0YAAIDf5OXlKSoqyufvazQE5ebmKjY2tty6hg0bKiwsTLm5ud5tWrZsWW6bsu/Jzc2tNASNHTtWo0aN8j4/ePCgEhMTtWPHDr/8IQaqgoICNW/eXD/++KMiIyNNl1Nt6Ju+nYC+6dsJ8vPz1aJFCzVq1Mgv73/WISgjI6PSUZaTZWVlqUOHDlV6v8oOZ1mWVW79z7cpGxar7Hslye12lxs5KhMVFeWoD0+ZyMhI+nYQ+nYW+nYWp/Zdq5Z/zuM66xCUnp6ugQMHnnabn4/cnEpcXJy++OKLcusOHDggj8fjHe2Ji4vzjgqV2bt3ryRVGEUCAACoqrMOQTExMYqJifHJD+/cubMmTZqknJwcxcfHS7InS7vdbrVv3967zR/+8AcVFxcrLCzMu01CQkKVwxYAAMDP+fU6QTt27NDq1au1Y8cOlZSUaPXq1Vq9erUOHz4sSerevbvatWun2267TatWrdLHH3+sMWPGaNiwYd7hvkGDBsntdmvIkCFat26d5s2bp8mTJ5/yzLDKuN1uPfnkk5UeIgtm9E3fTkDf9O0E9O2fvl2WH686OGTIEL3++usV1i9ZskRdu3aVZAel++67T5988onCw8M1aNAgTZ06tVzDa9eu1YgRI/Tll1+qYcOGuueee/TEE09UOQQBAAD8nF9DEAAAQKDi3mEAAMCRCEEAAMCRCEEAAMCRCEEAAMCRgioETZo0SV26dFFERIQaNGhQ6TY7duxQ3759VbduXcXExOiBBx5QcXFxuW3Wrl2rlJQUhYeHq2nTppowYYLfbt7mD0uXLpXL5ar0Kysry7tdZa+/8sorBiv/5Vq2bFmhp0cffbTcNlX5DNQk27Zt09ChQ5WUlKTw8HC1bt1aTz75ZIWegnF/z5gxQ0lJSapTp47at2+vzz77zHRJPvXUU0/p8ssvV/369dWkSRP169dPGzduLLfNkCFDKuzXTp06GarYNzIyMir0FBcX533dsixlZGQoISFB4eHh6tq1q9avX2+wYt+o7N8vl8ulESNGSAqeff3pp5+qb9++SkhIkMvl0nvvvVfu9ars36KiIt1///2KiYlR3bp1df3112vnzp1nXYvRe4f5WnFxsQYMGKDOnTtr5syZFV4vKSlR79691bhxYy1fvlx5eXkaPHiwLMvStGnTJNn3Z0lLS1NqaqqysrL0/fffa8iQIapbt65Gjx5d3S2dky5duignJ6fcuscff1yLFy+ucDuTzMxM9ezZ0/s8GO6tNmHCBA0bNsz7vF69et7HVfkM1DTfffedSktL9eqrr+q8887TunXrNGzYMB05ckRTp04tt20w7e+5c+dq5MiRmjFjhq688kq9+uqr6tWrlzZs2KAWLVqYLs8nli1bphEjRujyyy/X8ePHNW7cOHXv3l0bNmxQ3bp1vdv17NlTmZmZ3udlF5atyS688EItXrzY+zwkJMT7+JlnntHzzz+vWbNmqW3btpo4caLS0tK0ceNG1a9f30S5PpGVlaWSkhLv83Xr1iktLU0DBgzwrguGfX3kyBFdcskluuOOO3TTTTdVeL0q+3fkyJGaP3++5syZo+joaI0ePVp9+vRRdnZ2uc/KGVlBKDMz04qKiqqw/oMPPrBq1apl7dq1y7vu7bffttxut5Wfn29ZlmXNmDHDioqKso4dO+bd5qmnnrISEhKs0tJSv9fuD8XFxVaTJk2sCRMmlFsvyZo3b56ZovwkMTHReuGFF075elU+A8HgmWeesZKSksqtC7b9fcUVV1j33HNPuXXnn3++9eijjxqqyP/27t1rSbKWLVvmXTd48GDrhhtuMFeUHzz55JPWJZdcUulrpaWlVlxcnPX000971x07dsyKioqyXnnllWqqsHo8+OCDVuvWrb3/9wTjvv75v0tV2b8HDx60QkNDrTlz5ni32bVrl1WrVi3ro48+OqufH1SHw85k5cqVSk5OVkJCgnddjx49VFRUpOzsbO82KSkp5S7W2KNHD+3evVvbtm2r7pJ94v3339e+ffs0ZMiQCq+lp6crJiZGl19+uV555RWVlpZWf4E+NmXKFEVHR+vSSy/VpEmTyh0WqspnIBjk5+dXetflYNnfxcXFys7OVvfu3cut7969u1asWGGoKv/Lz8+XpAr7dunSpWrSpInatm2rYcOGee+vWJNt2rRJCQkJSkpK0sCBA/XDDz9IkrZu3arc3Nxy+97tdislJSWo9n1xcbHeeust3XnnneUuDByM+/pkVdm/2dnZ8ng85bZJSEhQcnLyWX8Ggupw2Jnk5uZWuOlqw4YNFRYW5r1Ja25uboV7kpV9T25urpKSkqqlVl+aOXOmevTooebNm5db/8c//lHdunVTeHi4Pv74Y40ePVr79u3TY489ZqjSX+7BBx/UZZddpoYNG+rLL7/U2LFjtXXrVv3lL3+RVLXPQE23ZcsWTZs2Tc8991y59cG0v/ft26eSkpIK+zI2NjZo9uPPWZalUaNG6Te/+Y2Sk5O963v16qUBAwYoMTFRW7du1eOPP65rrrlG2dnZNfYWCx07dtQbb7yhtm3bas+ePZo4caK6dOmi9evXe/dvZft++/btJsr1i/fee08HDx4s98trMO7rn6vK/s3NzVVYWJgaNmxYYZuz/fsf8CEoIyND48ePP+02WVlZFea6nEplt9qwLKvc+p9vY/00Kdr0bTrO5c9i586dWrBggd55550K2578n9+ll14qyZ5PE2j/KZ5N37///e+96y6++GI1bNhQN998s3d0SKraZyAQnMv+3r17t3r27KkBAwborrvuKrdtTdnfZ6Oyv6uBth99JT09XWvWrNHy5cvLrb/lllu8j5OTk9WhQwclJibqP//5j/r371/dZfpEr169vI8vuugide7cWa1bt9brr7/unQgc7Pt+5syZ6tWrV7lR62Dc16dyLvv3XD4DAR+C0tPTNXDgwNNuU9W7ycfFxemLL74ot+7AgQPyeDze1BkXF1chSZYNN/48mVa3c/mzyMzMVHR0tK6//vozvn+nTp1UUFCgPXv2GO/1ZL/kM1D2D+bmzZsVHR1dpc9AoDjbvnfv3q3U1FR17txZr7322hnfP1D3d1XExMQoJCSk0r+rNa2Xqrj//vv1/vvv69NPP1WzZs1Ou218fLwSExO1adOmaqrO/+rWrauLLrpImzZtUr9+/STZowHx8fHebYJp32/fvl2LFy/Wu+++e9rtgnFfl50FeLr9GxcXp+LiYh04cKDcaNDevXvVpUuXs/uB5zaVKbCdaWL07t27vevmzJlTYWJ0gwYNrKKiIu82Tz/9dI2cGF1aWmolJSVZo0ePrtL206ZNs+rUqVNuUnhNN3/+fEuStX37dsuyqvYZqIl27txptWnTxho4cKB1/PjxKn1PTd/fV1xxhXXvvfeWW3fBBRcE1cTo0tJSa8SIEVZCQoL1/fffV+l79u3bZ7ndbuv111/3c3XV59ixY1bTpk2t8ePHeyfOTpkyxft6UVFRUE2MfvLJJ624uDjL4/Gcdrtg2Nc6xcTo0+3fsonRc+fO9W6ze/fuc5oYHVQhaPv27daqVaus8ePHW/Xq1bNWrVplrVq1yjp06JBlWZZ1/PhxKzk52erWrZv19ddfW4sXL7aaNWtmpaene9/j4MGDVmxsrHXrrbdaa9eutd59910rMjLSmjp1qqm2ztnixYstSdaGDRsqvPb+++9br732mrV27Vpr8+bN1p///GcrMjLSeuCBBwxU6hsrVqywnn/+eWvVqlXWDz/8YM2dO9dKSEiwrr/+eu82VfkM1DS7du2yzjvvPOuaa66xdu7caeXk5Hi/ygTj/p4zZ44VGhpqzZw509qwYYM1cuRIq27duta2bdtMl+Yz9957rxUVFWUtXbq03H4tLCy0LMuyDh06ZI0ePdpasWKFtXXrVmvJkiVW586draZNm1oFBQWGqz93o0ePtpYuXWr98MMP1ueff2716dPHql+/vnffPv3001ZUVJT17rvvWmvXrrVuvfVWKz4+vkb3XKakpMRq0aKF9cgjj5RbH0z7+tChQ97/nyV5/90u+2W1Kvv3nnvusZo1a2YtXrzY+vrrr61rrrnGuuSSS6r8S2CZoApBgwcPtiRV+FqyZIl3m+3bt1u9e/e2wsPDrUaNGlnp6ekVfhNes2aNddVVV1lut9uKi4uzMjIyatwokGVZ1q233mp16dKl0tc+/PBD69JLL7Xq1atnRUREWMnJydaLL754xt88All2drbVsWNHKyoqyqpTp471q1/9ynryySetI0eOlNuuKp+BmiQzM7PSz/3JA73BuL8ty7L+9Kc/WYmJiVZYWJh12WWXlTt1PBicar9mZmZalmVZhYWFVvfu3a3GjRtboaGhVosWLazBgwdbO3bsMFv4L3TLLbdY8fHxVmhoqJWQkGD179/fWr9+vff10tJS72iJ2+22rr76amvt2rUGK/adBQsWWJKsjRs3llsfTPt6yZIllX6uBw8ebFlW1fbv0aNHrfT0dKtRo0ZWeHi41adPn3P6s3BZVg26FDIAAICPOOo6QQAAAGUIQQAAwJEIQQAAwJEIQQAAwJEIQQAAwJEIQQAAwJEIQQAAwJEIQQAAwJEIQQAAwJEIQQAAwJEIQQAAwJH+P1M3cdgLsAktAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import math\n",
    "from scipy.spatial.transform import Rotation as Rot\n",
    "import matplotlib.pyplot as plt\n",
    "def plot_covariance_ellipse(xEst, PEst):  # pragma: no cover\n",
    "    Pxy = PEst[0:2, 0:2]\n",
    "    eigval, eigvec = np.linalg.eig(Pxy)\n",
    "\n",
    "    if eigval[0] >= eigval[1]:\n",
    "        bigind = 0\n",
    "        smallind = 1\n",
    "    else:\n",
    "        bigind = 1\n",
    "        smallind = 0\n",
    "\n",
    "    t = np.arange(0, 2 * math.pi + 0.1, 0.1)\n",
    "    a = math.sqrt(eigval[bigind])\n",
    "    b = math.sqrt(eigval[smallind])\n",
    "    x = [a * math.cos(it) for it in t]\n",
    "    y = [b * math.sin(it) for it in t]\n",
    "    angle = math.atan2(eigvec[1, bigind], eigvec[0, bigind])\n",
    "    rot = Rot.from_euler('z', angle).as_matrix()[0:2, 0:2]\n",
    "    fx = rot @ (np.array([x, y]))*10\n",
    "    \n",
    "    print(xEst.shape)\n",
    "    \n",
    "    px = np.array(fx[0, :]+xEst[0,0]).flatten()\n",
    "    py = np.array(fx[1, :]+xEst[1,0]).flatten()\n",
    "    plt.scatter(xEst[0,0],xEst[1,0])\n",
    "    plt.plot(px, py, \"--r\")\n",
    "    plt.grid()\n",
    "    plt.xlim([-100,100])\n",
    "    plt.ylim([-100,100])\n",
    "    \n",
    "    \n",
    "xEst=np.ones([2,1])\n",
    "PEst=np.diag([5,40])\n",
    "plot_covariance_ellipse(xEst,PEst)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.64"
      ]
     },
     "execution_count": 43,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import numpy as np\n",
    "P=np.diag([4,1])\n",
    "\n",
    "S=np.diag([1,4])\n",
    "\n",
    "\n",
    "np.linalg.det(np.linalg.inv((np.linalg.inv(P)+np.linalg.inv(S))))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.6"
      ]
     },
     "execution_count": 44,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.trace(np.linalg.inv((np.linalg.inv(P)+np.linalg.inv(S))))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "6.45 µs ± 239 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)\n"
     ]
    }
   ],
   "source": [
    "%timeit np.linalg.det(np.random.randn(3,3))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "8.55 µs ± 178 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)\n"
     ]
    }
   ],
   "source": [
    "%timeit np.linalg.inv(np.random.randn(3,3))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {},
   "outputs": [],
   "source": [
    "A=np.random.randn(3,3)\n",
    "B=np.random.randn(3,3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[-1.22145477  1.05038512  1.71829976]\n",
      "[ 1.42659712  0.61873052 -1.61372014]\n",
      "[ 2.31043648 -0.3149285  -0.01667037]\n"
     ]
    }
   ],
   "source": [
    "print(np.linalg.eigvals(A))\n",
    "print(np.linalg.eigvals(B))\n",
    "print(np.linalg.eigvals(A+B))\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 109,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([9.  , 2.25])"
      ]
     },
     "execution_count": 109,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def rot(theta):\n",
    "    return np.array([[np.cos(theta),-np.sin(theta)],[np.sin(theta),np.cos(theta)]])\n",
    "\n",
    "P=np.diag([8,2])\n",
    "S=np.diag([4,1])\n",
    "\n",
    "t=np.deg2rad(60)\n",
    "\n",
    "np.linalg.eigvals(rot(t)*S*(rot(t).transpose())+P)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 110,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 1.  , -0.  ],\n",
       "       [-0.  ,  0.25]])"
      ]
     },
     "execution_count": 110,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "rot(t)*S*(rot(t).transpose())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 111,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.62281434, -0.00569803, -1.38849807],\n",
       "       [ 0.18208897,  0.7760501 ,  0.38214455],\n",
       "       [-0.24419824, -0.28542622, -0.68355794]])"
      ]
     },
     "execution_count": 111,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.linalg.inv(A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 112,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0.125, 0.   ],\n",
       "       [0.   , 0.5  ]])"
      ]
     },
     "execution_count": 112,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.linalg.inv(P)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 125,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "15.999999999999998"
      ]
     },
     "execution_count": 125,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "t=np.deg2rad(90)\n",
    "\n",
    "np.linalg.det(P+rot(t)*S*(rot(t).transpose()))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 126,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 1.49975978e-32, -0.00000000e+00],\n",
       "       [-0.00000000e+00,  3.74939946e-33]])"
      ]
     },
     "execution_count": 126,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "rot(t)*S*(rot(t).transpose())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 194,
   "metadata": {},
   "outputs": [],
   "source": [
    "from matplotlib.patches import Ellipse\n",
    "from matplotlib import pyplot as plt\n",
    "\n",
    "fig, ax = plt.subplots(1, 1, subplot_kw={\"aspect\": \"equal\"})\n",
    "\n",
    "angle = np.linspace(0, 135, 1)\n",
    "ellipse = [Ellipse((0, 0), 4, 2, a) for a in angle]\n",
    "\n",
    "for E in ellipse:\n",
    "    ax.add_patch(E)\n",
    "    E.set_edgecolor(\"black\")\n",
    "    E.set_facecolor(\"w\")\n",
    "    E.set_alpha(0.4)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 198,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.patches.Ellipse at 0x7ff0b69a1780>"
      ]
     },
     "execution_count": 198,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY0AAAD8CAYAAACLrvgBAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAGPVJREFUeJzt3X2MXXd95/H3156M7bEdPz+MHxI7i4lDVpsAo5QstEADwbHYBBa2NSu1WUByYUEq2l2psJEoAq1U2lJ221AiF6JCxfLQh5AIDI6hrABhHsbBCSFxEts4wfFzxg/x42Ts7/5xz5g743tnTjxz5s4475d0dc/D757z9e8en8/cc849NzITSZLKmNTqAiRJE4ehIUkqzdCQJJVmaEiSSjM0JEmlGRqSpNJGJTQi4t6IOBgRj9ZNmxsRmyPiqeJ5TpPX3lm0eSoi7hyNeiRJ1RitTxp/D6wZNO3DwHczcxXw3WJ8gIiYC/wp8FvATcCfNgsXSVLrjUpoZOb3gZ5Bk+8AvlAMfwF4W4OXvgXYnJk9mXkE2MzF4SNJGifaKlz2oszcB5CZ+yJiYYM2S4Ff143vKaZdJCLWA+sBpk+f/urVq1ePcrmSdHnbunXr4cxcMJJlVBkaZUSDaQ3va5KZG4ANAF1dXdnd3V1lXZJ02YmIp0e6jCqvnjoQEZ0AxfPBBm32AMvrxpcBeyusSZI0AlWGxgNA/9VQdwL3N2izCbg1IuYUJ8BvLaZJksah0brk9svAFuDaiNgTEe8F/gx4c0Q8Bby5GCciuiLicwCZ2QN8AvhZ8fh4MU2SNA7FRLw1uuc0JOnFi4itmdk1kmX4jXBJUmmGhiSpNENDklSaoSFJKs3QkCSVZmhIkkozNCRJpRkakqTSDA1JUmmGhiSpNENDklSaoSFJKs3QkCSVZmhIkkozNCRJpRkakqTSDA1JUmmGhiSptEpDIyKujYhtdY/jEfGhQW3eEBHH6tp8tMqaJEmXrq3KhWfmE8CNABExGXgWuK9B0x9k5lurrEWSNHJjeXjqFmBnZj49huuUJI2isQyNdcCXm8y7OSIejohvRcT1Y1iTJOlFGJPQiIh24HbgHxvMfgi4OjNvAP4G+HqTZayPiO6I6D506FB1xUqSmhqrTxq3AQ9l5oHBMzLzeGaeKIY3AldExPwG7TZkZldmdi1YsKD6iiVJFxmr0HgXTQ5NRcTiiIhi+KaipufGqC5J0otQ6dVTABHRAbwZ+KO6ae8DyMx7gHcC74+IPuA0sC4zs+q6JEkvXuWhkZmngHmDpt1TN3w3cHfVdUiSRs5vhEuSSjM0JEmlGRqSpNIMDUlSaYaGJKk0Q0OSVJqhIUkqzdCQJJVmaEiSSjM0JEmlGRqSpNIMDUlSaYaGJKk0Q0OSVJqhIUkqzdCQJJVmaEiSSjM0JEmlGRqSpNIqD42I2B0Rv4iIbRHR3WB+RMRfR8SOiHgkIl5VdU2SpEvTNkbreWNmHm4y7zZgVfH4LeCzxbMkaZwZD4en7gC+mDU/BmZHRGeri5IkXWwsQiOBByNia0SsbzB/KfDruvE9xbQBImJ9RHRHRPehQ4cqKlWSNJSxCI3XZuarqB2G+kBE/M6g+dHgNXnRhMwNmdmVmV0LFiyook5J0jAqD43M3Fs8HwTuA24a1GQPsLxufBmwt+q6JEkvXqWhERHTI2Jm/zBwK/DooGYPAH9YXEX1GuBYZu6rsi5J0qWp+uqpRcB9EdG/rv+bmd+OiPcBZOY9wEZgLbADOAW8u+KaJEmXqNLQyMxdwA0Npt9TN5zAB6qsQ5I0OsbDJbeSpAnC0JAklWZoSJJKMzQkSaUZGpKk0gwNSVJphoYkqTRDQ5JUmqEhSSrN0JAklWZoSJJKMzQkSaUZGpKk0gwNSVJphoYkqTRDQ5JUmqEhSSrN0JAklVZZaETE8oj4XkQ8HhG/jIg/btDmDRFxLCK2FY+PVlWPJGnkqvyN8D7gv2fmQxExE9gaEZsz87FB7X6QmW+tsA5J0iip7JNGZu7LzIeK4eeBx4GlVa1PklS9MTmnERErgFcCP2kw++aIeDgivhUR1w+xjPUR0R0R3YcOHaqoUknSUCoPjYiYAfwz8KHMPD5o9kPA1Zl5A/A3wNebLSczN2RmV2Z2LViwoLqCJUlNVRoaEXEFtcD4Umb+y+D5mXk8M08UwxuBKyJifpU1SZIuXZVXTwXweeDxzPyrJm0WF+2IiJuKep6rqiZJ0shUefXUa4E/AH4REduKaf8TuAogM+8B3gm8PyL6gNPAuszMCmuSJI1AZaGRmT8EYpg2dwN3V1WDJGl0+Y1wSVJphoYkqTRDQ5JUmqEhSSrN0JAklWZoSJJKMzQkSaUZGpKk0gwNSVJpVd5GRJqQMpNz585dePT19Q0YPn/+PJlJ/x1v6of7RQTFbdUujE+ePPnCo62t7aLh+vbSeGVo6LJy7tw5zp49e9Gjt7e36fS+vr4BwXDu3Lkhd+71O/jBz/3qA6X/+fz58xcFUP3wpEmTBqynra2N9vZ2pkyZwpQpUwYMD360t7czdepU2tr8L61quYVpwujt7eXkyZNDPnp7e5vuUNvb25k5cybz588fMO+KK64YEAyTJk1qyV/9gz/d9PX1XQi7+tA7efIkPT09DUNw8uTJdHR0MGPGDKZPn97wMXXq1DH/t+nyYWhoXDlz5gxHjx7l2LFjF56PHTvGiRMnyMyLdobz5s3jqquuGrBDnKiHefqDayT6Q6X+cejQIXbv3n1hvK+vj+nTp3PllVcye/ZsZs+ezaxZs5g9ezbTp08fpX+NLleGhlrm6NGjHDhwgJ6eHo4cOUJPTw99fX0DdmIve9nLmD17NjNmzKC9vb3VJY97/Z+g5s6d27RNX18fJ06c4Pjx4xw9epTDhw+zc+dOjhw5wvnz55k7d+6Fx/z585k/fz6TJnnNjGoMDY2JzKSnp4d9+/ZdeFxxxRUsWrSIefPmsXTpUubOncuMGTNaXeplr62t7cInjKuuumrAvNOnT18I8MOHD/PYY4/x/PPPs2jRIjo7O+ns7GTBggUj/kSkicvQUGV6e3vZsWMHzzzzDPv376ejo4POzk5WrlzJzTffbECMQ9OmTWPatGksWbLkwrSzZ8+yf/9+9u3bx5YtWzhy5AgLFixg2bJlrFq1yvfxJSYm4g/ldXV1ZXd3d6vLUBP79u1j+/btPP300yxfvpyVK1fS2dnJtGnTWl2aRsELL7zAgQMH2L17Nzt37mThwoWsXr2aq6++2sNY41xEbM3MrhEtw9DQaHnyySf5+c9/zqRJk1i9ejWrVq3ySp3LXF9fH7/61a/Yvn07R48e5frrr+fGG280PMapCREaEbEG+D/AZOBzmflng+ZPAb4IvBp4Dvj9zNw91DINjfHlhRde4Ac/+AE9PT287nWvY/Hixa0uSS1w9OhRtmzZQm9vL7fccouHrcah0QiNSv8ciIjJwGeA24BXAO+KiFcMavZe4Ehmvgz4NPDJKmvS6Nu0aROTJk3ibW97m4HxEjZ79mzWrFnDihUruP/+++nt7W11SapA1Z8hbwJ2ZOauzOwFvgLcMajNHcAXiuF/Am6JiXqh/UvQs88+y5kzZ3j961/vt5FFRHDDDTewdOlSHnvssVaXowpUHRpLgV/Xje8ppjVsk5l9wDFg3uAFRcT6iOiOiO5Dhw5VVK5erL6+Ptrb2yfsF+pUjfb2dvr6+lpdhipQdWg02pMMPolSpg2ZuSEzuzKza8GCBaNSnEZu+fLlnD17lm3btl100z69NO3evZsdO3Zw3XXXtboUVaDq0NgDLK8bXwbsbdYmItqAWUBPxXVplEyaNInbbruN3bt3s2nTJk6dOtXqktQiL7zwAlu2bOFHP/oRa9as8ZYkl6mqD0L/DFgVESuBZ4F1wH8e1OYB4E5gC/BO4F/TP1knlBkzZnD77bfT3d3N1772NZYsWcJ1113HsmXLPGz1EnDw4EG2b9/Orl27WL58Oe94xzuYMmVKq8tSRcbiktu1wP+mdsntvZn5vyLi40B3Zj4QEVOBfwBeSe0TxrrM3DXUMr3kdvzq7e1l586dbN++ndOnT3PttdeycuVK5syZY4BcRk6cOMHTTz/N448/Tl9fH6tXr+blL385HR0drS5NQ5gQ39OogqExMTz33HM88cQTPPPMM/T29rJ48eIL9y+aN2+eITKBHDt2bMB9w86dO8fSpUtZvXo1nZ2dvpcThKGhCePkyZMX7l+0b98+Tp48yeLFi1m0aBFz585lzpw5zJw5053POHDq1Cl6enro6enh4MGD7N+/n0mTJl0I/M7OTmbNmtXqMnUJDA1NWGfOnGH//v0Dbo1+9uxZ5syZM+DW6LNmzWLWrFneVXWUnT9/nueff37A75YcPXqUnp4eImLArdGXLFnit7svE6MRGn4bSy0xdepUVqxYwYoVKy5M6+3tpaen58JO7Mknn+TYsWMcP36cqVOnDvuLdH65sOb8+fPD/sLhqVOn6OjouHCL9Llz53LNNdcwZ84cz0toSP4v07jR3t7O4sWLL7oVSWZy+vRpTpw4MWDH19PTM2C8ra2Njo4Opk6d2vTnXhv9FOx4vbleZjb8udehHqdOneLs2bN0dHRcFKoLFy4cMD5e/90a3wwNjXsRQUdHx7B/AZ85c4ZTp05x5syZCzvR/p1ts9/U7u3tJSJoa2sb8Dvhww33n3tp9JyZF77oWP98/vx5+vr6BvwO+HDDzYKuf/jKK68cMK+jo4Np06Z5bkiVMTR02Zg6deol3Yq90c56uPHBwTD43GBEDNhxR8RFwVNm3J2/xhtDQy95/Ttof4NcGp4HNSVJpRkakqTSDA1JUmmGhiSpNENDklSaoSFJKs3QkCSVZmhIkkozNCRJpRkakqTSDA1JUmmV3HsqIv4C+A9AL7ATeHdmHm3QbjfwPHAO6Bvpj4NIkqpV1SeNzcC/zcx/BzwJfGSItm/MzBsNDEka/yoJjcx8MDP7itEfA8uqWI8kaWyNxTmN9wDfajIvgQcjYmtErB9qIRGxPiK6I6L70KFDo16kJGl4l3xOIyK+AyxuMOuuzLy/aHMX0Ad8qcliXpuZeyNiIbA5IrZn5vcbNczMDcAGgK6urmzURpJUrUsOjcx801DzI+JO4K3ALTn4Z81+s4y9xfPBiLgPuAloGBqSpNar5PBURKwB/gS4PTNPNWkzPSJm9g8DtwKPVlGPJGl0VHVO425gJrVDTtsi4h6AiFgSERuLNouAH0bEw8BPgW9m5rcrqkeSNAoq+Z5GZr6syfS9wNpieBdwQxXrlyRVw2+ES5JKMzQkSaUZGpKk0gwNSVJphoYkqTRDQ5JUmqEhSSrN0JAklWZoSJJKMzQkSaUZGpKk0gwNSVJphoYkqTRDQ5JUmqEhSSrN0JAklWZoSJJKMzQkSaVVFhoR8bGIeLb4jfBtEbG2Sbs1EfFEROyIiA9XVY8kaeQq+Y3wOp/OzL9sNjMiJgOfAd4M7AF+FhEPZOZjFdclSboErT48dROwIzN3ZWYv8BXgjhbXJElqourQ+GBEPBIR90bEnAbzlwK/rhvfU0y7SESsj4juiOg+dOhQFbVKkoYxotCIiO9ExKMNHncAnwX+DXAjsA/4VKNFNJiWjdaVmRsysyszuxYsWDCSsiVJl2hE5zQy801l2kXE3wHfaDBrD7C8bnwZsHckNUmSqlPl1VOddaNvBx5t0OxnwKqIWBkR7cA64IGqapIkjUyVV0/9eUTcSO1w027gjwAiYgnwucxcm5l9EfFBYBMwGbg3M39ZYU2SpBGoLDQy8w+aTN8LrK0b3whsrKoOSdLoafUlt5KkCcTQkCSVZmhIkkozNCRJpRkakqTSDA1JUmmGhiSpNENDklSaoSFJKs3QkCSVZmhIkkozNCRJpRkakqTSDA1JUmmGhiSpNENDklSaoSFJKs3QkCSVVsnPvUbEV4Fri9HZwNHMvLFBu93A88A5oC8zu6qoR5I0OioJjcz8/f7hiPgUcGyI5m/MzMNV1CFJGl2VhEa/iAjg94DfrXI9kqSxUfU5jd8GDmTmU03mJ/BgRGyNiPUV1yJJGqFL/qQREd8BFjeYdVdm3l8Mvwv48hCLeW1m7o2IhcDmiNiemd9vsr71wHqAq6666lLLliSNQGRmNQuOaAOeBV6dmXtKtP8YcCIz/3K4tl1dXdnd3T3yIiXpJSQito70gqMqD0+9CdjeLDAiYnpEzOwfBm4FHq2wHknSCFUZGusYdGgqIpZExMZidBHww4h4GPgp8M3M/HaF9UiSRqiyq6cy8780mLYXWFsM7wJuqGr9kqTR5zfCJUmlGRqSpNIMDUlSaYaGJKk0Q0OSVJqhIUkqzdCQJJVmaEiSSjM0JEmlGRqSpNIMDUlSaYaGJKk0Q0OSVJqhIUkqzdCQJJVmaEiSSjM0JEmlGRqSpNIMDUlSaSMKjYj4TxHxy4g4HxFdg+Z9JCJ2RMQTEfGWJq9fGRE/iYinIuKrEdE+knokSdUa6SeNR4H/CHy/fmJEvAJYB1wPrAH+NiImN3j9J4FPZ+Yq4Ajw3hHWI0mq0IhCIzMfz8wnGsy6A/hKZp7NzF8BO4Cb6htERAC/C/xTMekLwNtGUo8kqVptFS13KfDjuvE9xbR684Cjmdk3RJsLImI9sL4YPRsRj45SrVWaDxxudRHDmAg1gnWONuscXROlzmtHuoBhQyMivgMsbjDrrsy8v9nLGkzLS2jzmxmZG4ANRU3dmdnVrO14MRHqnAg1gnWONuscXROpzpEuY9jQyMw3XcJy9wDL68aXAXsHtTkMzI6ItuLTRqM2kqRxpKpLbh8A1kXElIhYCawCflrfIDMT+B7wzmLSnUCzTy6SpHFgpJfcvj0i9gA3A9+MiE0AmflL4GvAY8C3gQ9k5rniNRsjYkmxiD8B/ltE7KB2juPzJVe9YSR1j6GJUOdEqBGsc7RZ5+h6ydQZtT/4JUkant8IlySVZmhIkkobt6Ex0W5RUqxjW/HYHRHbmrTbHRG/KNqN+PK3S6jzYxHxbF2ta5u0W1P0746I+HAL6vyLiNgeEY9ExH0RMbtJu5b053D9U1wE8tVi/k8iYsVY1VZXw/KI+F5EPF78X/rjBm3eEBHH6raHj451nUUdQ76PUfPXRX8+EhGvGuP6rq3ro20RcTwiPjSoTcv6MiLujYiD9d9fi4i5EbG52Adujog5TV57Z9HmqYi4c9iVZea4fADXUfsiyv8DuuqmvwJ4GJgCrAR2ApMbvP5rwLpi+B7g/WNY+6eAjzaZtxuY38J+/RjwP4ZpM7no12uA9qK/XzHGdd4KtBXDnwQ+OV76s0z/AP8VuKcYXgd8tQXvdSfwqmJ4JvBkgzrfAHxjrGt7se8jsBb4FrXvd70G+EkLa50M7AeuHi99CfwO8Crg0bppfw58uBj+cKP/Q8BcYFfxPKcYnjPUusbtJ42coLcoKdb9e8CXx2J9FbkJ2JGZuzKzF/gKtX4fM5n5YP7mbgE/pvY9nvGiTP/cQW27g9p2eEuxbYyZzNyXmQ8Vw88DjzPEXRfGuTuAL2bNj6l9x6uzRbXcAuzMzKdbtP6LZOb3gZ5Bk+u3wWb7wLcAmzOzJzOPAJup3S+wqXEbGkNYCvy6bnzEtygZZb8NHMjMp5rMT+DBiNha3BqlFT5YfMS/t8lH1jJ9PJbeQ+2vzEZa0Z9l+udCm2I7PEZtu2yJ4vDYK4GfNJh9c0Q8HBHfiojrx7Sw3xjufRxP2+Q6mv9ROB76st+izNwHtT8ggIUN2rzofq3q3lOlxDi5RUlZJet9F0N/ynhtZu6NiIXA5ojYXvyVMGqGqhP4LPAJav3xCWqH0t4zeBENXjvq12aX6c+IuAvoA77UZDGV92cDLdsGL0VEzAD+GfhQZh4fNPshaodZThTnt75O7cu4Y22493Fc9GdxbvR24CMNZo+XvnwxXnS/tjQ0coLdomS4eiOijdqt4l89xDL2Fs8HI+I+aoc6RnUnV7ZfI+LvgG80mFWmj0esRH/eCbwVuCWLA7ANllF5fzZQpn/62+wptotZXHz4oHIRcQW1wPhSZv7L4Pn1IZKZGyPibyNifmaO6c33SryPY7JNlnAb8FBmHhg8Y7z0ZZ0DEdGZmfuKQ3kHG7TZQ+1cTL9l1M4jNzURD0+N51uUvAnYnpl7Gs2MiOkRMbN/mNrJ3jG9W++g48Bvb7L+nwGronYFWju1j+MPjEV9/SJiDbU7BtyemaeatGlVf5bpnweobXdQ2w7/tVnwVaU4h/J54PHM/KsmbRb3n2uJiJuo7ROeG7sqS7+PDwB/WFxF9RrgWP+hlzHW9EjCeOjLQeq3wWb7wE3ArRExpzhUfWsxrblWnOkveTXA26ml4FngALCpbt5d1K5eeQK4rW76RmBJMXwNtTDZAfwjMGUMav574H2Dpi0BNtbV9HDx+CW1wzBj3a//APwCeKTYqDoH11mMr6V2tc3OFtW5g9qx1m3F457BdbayPxv1D/BxaiEHMLXY7nYU2+E1LejD11E71PBIXT+uBd7Xv50CHyz67mFqFxz8+xbU2fB9HFRnAJ8p+vsX1F1ROYZ1dlALgVl108ZFX1ILsn3AC8V+873UzqF9F3iqeJ5btO0CPlf32vcU2+kO4N3DrcvbiEiSSpuIh6ckSS1iaEiSSjM0JEmlGRqSpNIMDUlSaYaGJKk0Q0OSVNr/B5FOipH7R+7zAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots()\n",
    "ax.axis([-10, 10, -10, 10])\n",
    "\n",
    "E=Ellipse((0,0),8,2,0)\n",
    "E.set_edgecolor(\"black\")\n",
    "E.set_facecolor(\"w\")\n",
    "E.set_alpha(0.4)\n",
    "ax.add_patch(E)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
